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NOMENCLATURE 

A  =  column  cross-sectional  area 

a  =  exponential  model  coefficients 

B  =  Bayes  factor 

b°  =  body  force  vector  with  respect  to  the  undefonned  configuration,  physical  coordinates 
b  =  buckling  state 

C  =  fourth-order  elasticity  tensor,  physical  coordinates 

C  =  Bayesian  hypothesis  testing-based  confidence  metric 
Cp  =  coefficient  of  pressure 

c  =  Eckert’s  reference  temperature  discrepancy  model  coefficients 

D  =  damping  matrix,  modal  coordinates 
D  =  diameter  of  the  spherical  dome  wind  tunnel  specimen 
d  =  heat  transfer  discrepancy  model  coefficients 

E  =  Green  strain  tensor,  physical  coordinates 

E  =  modulus  of  elasticity 

e  =  model  error 

F  =  deformation  gradient  tensor,  physical  coordinates 

Fcr  =  Euler  buckling  load  for  a  column 
Fj  =  ith  force  from  external  excitation,  modal  coordinates 
f\  =  first  natural  frequency 

fs  =  third  natural  frequency 

FI  =  height  of  the  spherical  dome  wind  tunnel  specimen 

/  =  moment  of  inertia 

Kn'  =  linear  stiffness  matrix,  modal  coordinates 

Kff  =  element  of  the  third-order  tensor  associated  with  quadratic  terms,  modal  coordinates 

=  element  of  the  fourth-order  tensor  associated  with  the  cubic  terms,  modal  coordinates 

K  =  number  of  possible  states 

L  =  column  length 

M  =  mass  matrix,  modal  coordinates 

M  =  Mach  number 

N  =  number  of  samples 

p  =  aerodynamic  pressure 

Q  =  aerodynamic  heat  flux 

q  =  dynamic  pressure  (pU  / 2) 

Req  =  equivalence  ratio  (fuel-to-air  ratio  divided  by  stoichiometric  fuel-to-air  ratio) 

S  =  second  Piola-Kirchhoff  stress  tensor,  physical  coordinates 

S  =  main  effect  sensitivity  index 

St  =  total  effect  sensitivity  index 

T  =  temperature 

Tbuck  =  critical  buckling  temperature 

t  =  time 

lfn>  =  modal  basis  function 

U  =  flow  velocity 

u  =  displacement  vector,  physical  coordinates 

W\  =  transverse  component  of  mode  shape  i 
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w 

= 

transverse  panel  displacement 

X 

= 

model  prediction,  location  along  dome 

y 

= 

observed  data 

c 

Z 

= 

critical  state  boundary  value 

z 

= 

vector  of  measured  variables 

a 

= 

coefficient  of  thennal  expansion 

P 

= 

oblique  shock  angle  relative  to  freestream 

S 

= 

tolerance  limit 

Sij 

= 

Kronecker  delta 

8 

= 

model  error,  measurement  error 

r 

= 

ratio  of  specific  heats 

* 

= 

uncertain  input  and  model  parameters 

P 

= 

mean 

n(0 

= 

response  vector,  modal  coordinates 

n 

= 

model  prediction 

Qo 

= 

structure  domain  in  the  undeformed  configuration 

n, 

= 

set  of  parents  for  node  i  in  a  Bayesian  network 

n 

= 

probability  density  function 

e 

= 

panel  inclination  angle  to  freestream 

p 

= 

density 

a 

= 

standard  deviation 

Subscripts 

aw 

= 

adiabatic  wall 

e 

= 

edge  of  boundary  layer 

i 

= 

variable  index,  position  along  spherical  dome 

true 

= 

true  value  of  x 

pred 

model  prediction  of  x 

w 

= 

wall,  aerodynamic  surface 

ERT 

= 

Eckert’s  reference  temperature  method 

OS 

= 

oblique  shock  relations 

PI 

= 

Piston  theory 

1 

freestream  flow 

3 

= 

flow  at  the  leading  edge  of  panel 

4 

= 

flow  at  location  of  interest  along  the  panel 

Superscripts 

fP 

= 

flat  plate 

sd 

= 

spherical  dome 

* 

= 

flow  properties  evaluated  at  Eckert’s  reference  temperature 
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ABSTRACT 

Lack  of  confidence  in  structural  response  and  life  predictions  of  a  vehicle  exposed  to 
combined  extreme  environments  has  consistently  prevented  the  USAF  from  fielding  affordable, 
reliable,  and  reusable  hypersonic  space  access  platforms.  Significant  strides  have  been  made  in 
modeling  complex  interactions  of  the  multi-physics,  fluid-thermal-structural  coupling  applicable 
to  hypersonic  flow  conditions.  However,  validation  of  these  models  remains  a  challenge  due  to 
limited  experimental  data  for  hypersonic  conditions.  This  research  addresses  fundamental  and 
critical  issues  in  quantifying  uncertainty  and  assessing  the  confidence  in  model  predictions  of 
hypersonic  structural  response  through  a  systematic  framework.  The  first  year  of  this  research 
focused  on  identifying  and  developing  the  components  of  the  model  uncertainty  framework  for 
aerodynamic  pressure  and  heating  predictions,  including  global  sensitivity  analysis,  Bayesian 
model  calibration,  and  validation  metric  comparison.  The  second  year  emphasized  effectively 
integrating  information  into  the  coupled  system  through  segmented  Bayesian  model  calibration. 
The  final  year  brought  together  model  discrepancy  in  aerodynamic  pressure  calibrated  from 
aerothermal  experiments  with  nonlinear  structural  dynamic  reduced  order  models  to  investigate 
the  uncertainty  and  sensitivity  in  coupled  aeroelastic  response  (i.e.,  flutter  and  limit  cycle 
oscillation). 
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1  INTRODUCTION 

Uncertainty  inherently  exists  in  all  computational  model  predictions  due  to  imperfect 
knowledge  and  physical  variability  in  the  system,  model  order  reduction,  assumptions  and 
approximations,  and  the  limited  experimental  data  available  for  model  validation.  This  is 
especially  the  case  for  structures  in  hypersonic  environments  due  to  the  complex  and  poorly- 
understood  loading  from  the  inherently  coupled  multi-physics  nature  of  the  fluid-thermal- 
structural  interaction.  In  traditional  detenninistic  design,  a  margin  of  safety  is  introduced  to 
safeguard  against  uncertainty.  However,  this  can  lead  to  an  inefficient  or  an  unrealizable  design. 
For  an  aircraft  to  achieve  the  demanding  perfonnance  requirements  of  sustained  hypersonic 
flight,  weight  is  a  very  significant  design  constraint  [1-3].  Therefore,  this  research  effort  was 
focused  on  acquiring  the  fundamental  understanding  required  for  integrating  various  sources  of 
uncertainty  in  a  coupled  aerothermoelastic  simulation,  identifying  the  most  significant  error 
sources,  and  developing  a  method  for  dynamically  quantifying  model  prediction  confidence 
during  transient,  combined,  aerothermal  and  aero-pressure  loading. 

Substantial  research  has  been  performed  on  investigating  the  model  components  for  the 
physics  of  a  coupled  aerothennoelastic  panel  and  the  solution  procedures  for  both  quasi-static 
and  dynamic  solutions  [4-12].  However,  the  current  state  of  the  art  focuses  on  detenninistic 
calculations  with  limited  uncertainty  analysis.  Lamorte  et  al.  investigated  the  implementation  of 
a  stochastic  collocation  approach  for  propagating  uncertainty  in  aerothennoelastic  analysis  [13]. 
Related  work  expanded  on  uncertainty  propagation  in  aerothermoelastic  analysis  for  hypersonic 
vehicles  with  emphasis  on  assessing  the  impact  of  aerothennoelastic  defonnation  on 
aerodynamic  heating  [14],  Culler  et  al.  also  identified  2-way  coupling  between  structural 
defonnation  and  aerodynamic  heating  as  an  important  consideration  in  modeling  an 
aerothermoelastic  panel  [11].  Ostoich  et  al.  looked  at  the  heat  flux  into  a  spherical  dome 
protuberance  on  a  flat  plate  model,  calculated  from  high-fidelity,  fully  compressible  Navier- 
Stokes  equations  without  turbulence  model  and  compared  the  results  to  experimental  data  and 
lower-order  methods  [15,16].  Rangavajhala  et  al.  investigated  the  discretization  error  associated 
with  multidisciplinary  analyses  caused  by  mesh  sizes  and  mismatch  of  disciplinary  meshes  [17]. 
These  efforts  underscore  the  importance  of  understanding  the  uncertainty  in  a  coupled 
aerothermoelastic  model;  however,  many  questions  remain  about  the  significant  sources  of 
uncertainty  and  how  to  assess  the  confidence  in  multi-physics  model  predictions. 

The  model  uncertainty  framework  used  in  this  research  is  founded  in  Bayesian  statistics, 
which  is  an  effective  approach  for  uncertainty  reduction  and  confidence  assessment  through  the 
integration  of  computational  prediction  and  experimental  observation.  For  highly  coupled,  multi¬ 
physics  problems  in  which  data  can  be  sparse  and  models  can  be  numerous,  it  is  necessary  to 
connect  them  in  a  systematic  way  to  reduce  uncertainty.  Bayesian  networks  are  used  to  reflect 
the  complex  relationships  between  sources  of  uncertainty  and  model  predictions,  which  are 
represented  as  nodes  in  the  network.  The  value  of  Bayesian  networks  lies  in  their  ability  to  apply 
experimental  data  to  individual  nodes  and  reduce  uncertainty  over  the  entire  network  [18,19]. 
This  is  particularly  useful  for  coupled,  aerothermoelastic  models  in  which  data  is  not  necessarily 
available  to  validate  the  fully-coupled  prediction;  however,  data  from  a  subset  of  the  coupled 
physics  can  be  readily  integrated  into  the  network. 

The  next  subsection  discusses  the  developed  model  uncertainty  framework  and  the  and  the 
aspects  of  uncertainty  quantification  (UQ)  and  validation  and  verification  (V&V)  considered 
during  this  research  effort. 
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1.1  Uncertainty  Quantification  and  Model  Prediction  Confidence  Assessment 

Quantifying  the  confidence  in  model  prediction  consists  of  two  intertwined,  yet  distinct, 
activities:  uncertainty  quantification  (UQ)  and  verification  and  validation  (V&V).  The  science  of 
uncertainty  quantification  for  numerical  simulations  (i.e.,  the  quantitative  characterization  and 
reduction  of  uncertainties)  has  origins  dating  back  to  the  early- 1 990 ’s.  However,  over  the  past 
decade  there  has  been  a  surge  in  multiple  research  communities  towards  fonnalizing  and 
generalizing  the  process.  Thus,  there  are  now  multiple  descriptions  and  implementations  of  the 
UQ  and  V&V  processes  that  are  all  similar  in  their  objectives,  but  different  in  their  details.  Due 
to  the  nature  of  UQ  and  V&V  research,  it  is  necessary  to  establish  terminology  and  research 
scope.  Uncertainty  is  inherent  in  all  computational  model  predictions  due  to  imperfect 
knowledge  and  physical  variability  in  the  system,  model  order  reduction,  assumptions  and 
approximations,  and  the  limited  experimental  data  available  for  model  calibration  and  validation. 
This  is  especially  the  case  for  compliant  structures  in  hypersonic  environments  due  to  the 
complex  and  poorly-understood  loading  and  the  coupled  multi-physics  nature  of  the  fluid- 
thermal-structural  interaction.  Physical  variability  is  incorporated  in  fluid-thennal-structural 
models  through  variations  in  material  properties,  geometry,  boundary  conditions,  and  load 
interactions.  The  aerothennoelastic  model  prediction  also  has  both  model-form  error  and 
numerical  errors.  Model-form  error  encompasses  the  errors  in  representing  the  physical  system 
with  a  particular  model.  These  errors  are  assessed  by  comparing  model  predictions  to  physical 
observations.  Numerical  errors  include  errors  from  sampling,  discretization,  coupled  solution 
procedures,  and  solution  approximation  error  from  model  order  reduction.  This  research  focused 
on  uncovering  the  most  significant  contributors  to  the  overall  uncertainty  from  each  individual 
model  component  of  the  coupled  system,  in  addition  to  quantifying  the  uncertainty  associated 
with  the  degree  of  coupling. 

Regarding  V&V,  the  primary  interest  for  this  research  is  on  validation,  due  to  the  significant 
challenges  posed  by  the  limited  experimental  data  available  for  structures  operating  in 
hypersonic  environments.  Validation  metrics  can  be  used  in  several  different  ways  for 
uncertainty  management,  including  model  selection,  model  validation,  or  model  prediction 
confidence.  Note  that  in  this  context,  model  prediction  confidence  refers  to  the  case  when  a 
validation  metric  must  be  extrapolated  beyond  the  validation  domain  of  the  experiment. 

Considering  the  UQ  and  V&V  discussion  above,  a  framework  for  quantifying  model 
prediction  confidence  for  hypersonic  aerostructures  was  developed  and  is  shown  in  Figure  1.1. 
The  first  component  focuses  on  characterizing  model  uncertainty,  which  is  accomplished  by 
constructing  a  Bayesian  network  and  calibrating  the  uncertain  parameters  and  model  discrepancy 
[20-22].  This  is  not  just  about  quantifying  the  uncertainty  in  each  of  the  individual  uni- 
disciplinary  constituent  models,  but  rather  about  quantifying  the  uncertainty  in  the  coupled 
interactions  of  the  multi-fidelity,  multi-physics  models,  as  well.  The  second  component 
addresses  efficiently  propagating  the  model  uncertainties  forward  to  a  quantity  of  interest  (Qol) 
in  the  coupled  system.  In  this  construct,  the  sensitivity  of  the  Qol  to  each  of  the  uncertain  inputs 
and  model  errors  can  be  evaluated,  which  is  important  for  identifying  where  to  focus  the 
resources  for  uncertainty  reduction  efforts.  The  third  and  final  component  is  managing  the 
uncertainty  by  using  data  for  multiple  different  purposes,  such  as  prediction  confidence 
assessment,  model  validation,  and/or  model  selection. 
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Figure  1.1.  Framework  for  integrating  of  models,  data,  and  uncertainty 


1.2  Report  Overview  and  Outline 

This  Lab  Task  effort  takes  the  initial  steps  in  a  longer-term  research  plan  for  developing  a 
framework  for  integrating  various  sources  of  uncertainty  in  a  coupled  hypersonic  structural 
simulation  and  assessing  the  confidence  in  model  predictions.  The  FY12-14  Lab  Task  consisted 
of  three  primary  objectives.  The  first  objective  is  to  develop  a  systematic  framework  for 
enabling  the  integration  of  the  various  sources  of  uncertainty  from  the  individual  disciplines  of  a 
fully-coupled  system.  The  second  objective  is  to  investigate  uncertainty  quantification  and 
analysis  of  fluid-thermal-structural  interactions  in  hypersonic  flow.  The  third  objective  is  to 
provide  a  decision-making  metric  for  determining  the  necessary  model  fidelity  and  degree  of 
coupling  to  achieve  the  desired  level  of  confidence  in  the  system  prediction.  All  three  objectives 
play  a  critical  role  in  assessing  the  predictive  capability  of  a  coupled  aerothennoelastic  system 
model.  The  objectives  were  addressed  through  several  related  investigations,  presented  in 
Sections  3-8  of  this  final  report. 

Section  2  introduces  the  aerothennoelastic  coupled  system  considered  throughout  this 
research,  including  aerothermal  models  and  historic  high-speed  wind  test  data.  The  first  year  of 
this  research  (Section  3)  focused  on  identifying  and  developing  the  components  of  the  model 
uncertainty  framework  for  aerodynamic  pressure  and  heating  predictions,  including  global 
sensitivity  analysis,  Bayesian  model  calibration,  and  confidence  assessment  [23,31].  The  second 
year  emphasized  effectively  integrating  information  into  the  coupled  system  through  segmented 
Bayesian  model  calibration,  as  well  as  using  time-dependent  aerothermal  data  in  the  Bayesian 
network,  discussed  in  Sections  4  and  5,  respectively  [24,25,33].  A  pre-validation  study  was  also 
conducted  (Section  7)  for  the  historic  aerothermal  data  to  determine  if  the  models  being 
considered  adequately  captured  the  observed  flow  characteristics  [27].  The  third  year  brought 
together  model  discrepancy  in  aerodynamic  pressure  calibrated  from  aerothermal  experiments 
with  nonlinear  structural  dynamic  reduced  order  models  to  investigate  the  uncertainty  and 
sensitivity  in  coupled  aeroelastic  response  (i.e.,  flutter  and  limit  cycle  oscillation),  which  is 
covered  in  Section  6  [26,32,34],  Finally,  many  real  world  systems  have  discrete  operating  states, 
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e.g.,  pre-buckled  or  post-buckled  columns,  laminar  or  turbulent  flow,  elastic  or  plastic  response, 
and  various  stages  of  wear.  Identification  of  precise  boundaries  between  states  is  typically  very 
challenging  due  to  incomplete  understanding  of  physical  processes  and  experimental  uncertainty. 
In  Section  8,  a  method  was  developed  for  quantifying  the  uncertainty  in  state  boundaries  for  thin 
beam  buckling  [28,35,36]. 
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2  AEROTHERMOELASTIC  COUPLING 

Aircraft  structures  exposed  to  extreme  environments  are  subjected  to  coupled  aerodynamic, 
thermal,  and  acoustic  loading  [4-12].  Neglecting  these  interactions  can  lead  to  gross  errors  in 
model  predictions  [10-13, 38, 39]. To  account  for  these  uncertainties,  hypersonic  aircraft  structures 
must  be  designed  to  within  accurately-quantified  safety  margins  that  are  not  overly  conservative, 
so  that  the  platform  can  have  the  minimum  structural  weight  that  pennits  proper  execution  of  the 
mission  objectives  [40],  When  the  structure  is  designed  to  these  margins,  it  must  necessarily 
operate  at  the  intersection  of  the  applicable  technical  disciplines  associated  with  extreme 
hypersonic  environments,  and  be  able  to  withstand  intense,  coupled,  structural,  fluid,  thennal, 
and  acoustic  loads  [4-8], 

Consider  a  panel  section  on  the  forebody  of  a  representative  hypersonic  vehicle  configuration, 
as  shown  in  Figure  2.1  [11].  As  the  vehicle  is  subjected  to  a  hypersonic  flow,  an  attached  oblique 
shock  is  created  at  the  forebody  leading  edge  (location  ‘1’).  This  results  in  aerodynamic  pressure 
applied  to  the  area  of  interest  (location  ‘4’),  causing  elastic  deformation  of  the  panel  into  the 
flow  field,  which  feeds  back  to  affect  the  aerodynamic  loads  on  the  panel.  This  is  commonly 
referred  to  as  the  aeroelastic  portion  of  the  coupling.  The  panel  is  also  subjected  to  aero  thennal 
effects  from  aerodynamic  heating.  This  aerothermal  component  is  also  coupled  to  the  aeroelastic 
component,  since  a  change  in  the  temperature  of  the  structure  causes  additional  deformation, 
which  in  turn  further  affects  both  the  aerodynamic  pressure  and  the  aerodynamic  heating.  The 
wetted  surface  of  the  structure  acts  as  a  boundary  condition  to  the  flow  problem  and  the 
aerodynamics  load  the  structure.  The  panel  experiences  extreme  aerodynamic  pressure  and 
heating  as  the  temperature  of  the  air  increases  within  the  inviscid  boundary  layer.  Heat  transfers 
to  the  body,  where  the  temperature  gradient  in  the  panel  induces  deformation  due  to  temperature- 
dependent  material  properties  and  thennal  moments,  in  addition  to  structural  dynamic 
deformation. 


Figure  2.1.  Representative  hypersonic  vehicle  structure  with  aerothermoelastic  panel  [11] 
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Figure  2.1  schematically  illustrates  the  fluid-thermal-structural  interactions  as  a  coupled 
aerothennoelastic  response,  consisting  of  the  following  individual  model  components: 
aerodynamic  pressure,  aerodynamic  heating,  heat  transfer,  and  structural  deformation.  Modeling 
these  coupled  disciplines  is  crucial  to  accurately  predicting  the  structural  response  under 
hypersonic  flow  conditions.  McNamara  and  Friedmann  (2007)  provide  a  comprehensive  review 
of  the  current  state-of-the-art  for  solution  strategies  for  calculating  the  response  of  a  hot  structure 
in  a  hypersonic  flow  [4]. 


Figure  2.2.  Coupling  in  aerothermoelasticity 

Up  to  this  point,  the  coupling  of  aero-thermo-structural  disciplines  has  been  discussed 
generally.  However,  the  fidelity  at  which  each  of  the  individual  disciplines  is  modeled,  and  how 
these  multi-fidelity  models  are  then  coupled  together,  is  a  significant  consideration  for  correctly 
quantifying  the  uncertainty  in  the  coupled  system.  The  aerothennoelastic  plate  can  be  modeled  at 
multiple  levels  of  fidelity  for  each  of  the  aerodynamic,  structural,  and  thermal  effects, 
respectively.  Limitations  on  computational  resources  and  the  acceptable  amount  of  uncertainty  in 
the  quantity  of  interest  are  often  competing  constraints  that  drive  the  degree  of  model  fidelity  and 
the  level  of  coupling  selected  for  simulation  purposes.  For  example,  the  aerodynamic  pressure 
and  heating  could  be  calculated  using  a  computational  fluid  dynamics  (CFD)  code  (e.g., 
CFL3D).  Conversely,  the  decision  could  be  made  to  use  lower-fidelity  aerothermal  model 
components,  such  as  piston  theory  and  Eckert’s  reference  methods  [9,1 1,15,12,41,42].  For  quasi¬ 
static  calculations,  the  CFD  approach  adds  minimal  cost.  However,  the  computational  time 
drastically  increases  for  a  transient  dynamic  solution.  To  reduce  the  computational  cost,  the 
aeroelastic  components  of  the  coupled  problem  can  be  approximated  with  nonlinear  reduced 
order  models  (ROMs)  built  from  finite  element  models  and  piston  theory  [43].  However,  the  use 
of  reduced  order  models  introduces  solution  approximation  and  model-form  errors.  Therefore,  in 
many  cases,  to  quantify  the  uncertainty  of  the  coupled  aerothermoelastic  system,  the  numerical 
and  model-form  errors  of  each  model  component  must  be  quantified  and  propagated  throughout 
the  system.  For  example,  model  uncertainty  in  aero-pressure  predictions  from  piston  theory 
could  be  the  primary  contributor  to  the  uncertainty  in  aerodynamic  heating  [24]  and  structural 
response  predictions  [26]  when  coupled  together.  These  are  two  of  the  key  observations  from 
this  Lab  Task  investigation,  which  will  be  covered  in  Sections  4  and  6,  respectively. 
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The  next  section  defines  several  of  the  model  components  in  the  aerothermoelastic  system 
shown  in  Figure  2.2. 

2.1  Aerothermoelastic  Models 

Given  the  freestream  flight  conditions  (p i,  M\,  T\)  and  the  surface  inclination  angle  ( 9 )  from 
Figure  2.1,  the  local  conditions  at  the  leading  edge  of  the  panel  (pi,  Mi,  Ti )  resulting  from  an 
oblique  shockwave  can  be  computed  using  oblique  shock  relations,  shown  in  Eqs.  (2.1)  -  (2.4). 
The  oblique  shock  calculations  do  not  have  any  dependency  on  the  geometry  of  the  panel  itself, 
solely  the  surface  inclination  of  the  forebody  and  freestream  conditions.  Thus,  it  is  only  valid  at 
locations  where  no  structural  defonnation  is  present  (i.e.  flat  plate). 


Pl  =1+  ly  (Mf sin2/?  1) 

P\  r+ 1 

(2.1) 

Pi _  (j  +  \)M  ,2  sin2(/?) 

A  (y-l)M12sin2(y0)  +  2 

(2.2) 

T3  _p3/pl 

T\  Pi  /  A 

M2  sin2  (/?)  H — 

M;  sin2  (ft  -  0)  =  — - 

(— i)M12sin2(y0)-l 

Y~  1 

(2.3) 

(2.4) 

Once  the  flow  properties  at  the  leading  edge  of  the  panel  are  calculated  from  oblique  shock 
relations,  piston  theory  provides  a  simplified  relationship  between  the  unsteady  pressure  on  the 
panel  and  turbulent  surface  pressure  [41].  This  simple  model  is  desired  for  computational 
tractability  and  uses  the  leading  edge  conditions  to  approximate  the  aerodynamic  pressure  load 
chord-wise  across  the  panel  (p4,  M4,  T4).  In  piston  theory,  the  pressure  prediction  is  dependent  on 
the  slope  of  the  panel  (dw/dx)  and  the  velocity  of  defonnation  (ow/ot).  A  3ld-order  expansion  of 
piston  theory  is  presented  in  Eq.  (2.5). 


.  q3  1  dw  dw 

p4=Pi  +  2-^- - +  —  + 

M3  ^  U3  dt  dx  j 


, ,  1  dw  dw 

^  U3  dt  dx 


J_du>+3wV  (25) 
U3  dt  dx 


After  calculating  the  aerodynamic  pressure  and  flow  conditions  along  the  panel  surface,  the 
aerodynamic  heat  flux  is  predicted  using  the  computationally  efficient  Eckert's  reference 
temperature  method  assuming  a  calorically  perfect  gas  [42].  The  Eckert's  reference  temperature 
is  computed  by  Eq.  (2.6)  and  the  heat  flux  across  the  spherical  dome  follows  in  Eq.  (2.7). 

r  =  T3  +  0.5(7;  -Te)  +  0.22(Taw  -  T3)  (2.6) 
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(2.7) 


Qa  =  St*  p*U ecp(Taw-Tw) 

where,  St*  is  the  reference  Stanton  number,  p  is  the  reference  density,  Ue  is  the  inviscid  flow 
velocity  at  the  dome  location,  c*  is  the  reference  specific  heat,  Taw  and  Tw  are  the  adiabatic  wall 
and  actual  wall  temperature,  respectively  and  Te  is  the  boundary  layer  edge  temperature. 

2.2  Aerothermal  Experiments 

Accurately  modeling  aero-thermo-elastic  response  for  hypersonic  aircraft  structures  is 
obviously  challenging,  especially  due  to  the  inability  to  replicate  the  bulk  of  the  in-flight  loading 
conditions  through  ground  test  facilities  and  hence  improve  the  models  being  used.  Laboratory 
experiments  to  approximate  these  unfamiliar  flight  environments  is  the  other  half  of  ensuring 
accurate  aero-thermo-acoustic  modeling.  However,  the  available  experimental  techniques  and 
facilities  are  not  fully  capable  of  simultaneously  capturing  coupled  aerothennoelastic  response  in 
hypersonic  flow.  Therefore,  we  must  rely  on  maximizing  the  utility  of  whatever  data  can  be 
acquired,  or  may  already  be  available,  for  a  subset  of  the  multi-physics  interactions.  For 
example,  historical  tests  performed  by  Glass  and  Hunt  in  1986  at  NASA’s  8ft  High-Temperature 
Wind  Tunnel  (HTT)  investigated  the  thennal  and  structural  loads  on  body  panels  in  hypersonic 
environments  [44].  These  tests  measured  the  aerodynamic  pressure  and  heating  on  spherical 
dome  protuberances  into  the  flow  that  simulated  deformed  aircraft  panels.  But  these  tests  were 
performed  on  rigid  domes  protruding  into  the  flow  and  were  not  instrumented  to  measure 
structural  response,  thus  the  dynamic  structural  response  was  not  captured.  The  8-foot  High- 
Temperature  Tunnel  can  simulate  up  to  Mach  7  flow  at  an  altitude  between  25  and  40  km  for  up 
to  2  minutes  by  combusting  a  mixture  of  methane  and  air.  The  flow  conditions  for  the  tests  of 
interest  had  a  turbulent  boundary-layer  at  the  panel  location,  and  the  panel  holder  had  a  sharp 
leading  edge,  similar  to  the  representative  hypersonic  vehicle  depicted  in  Figure  2.1. 

The  experiments  performed  by  Glass  and  Hunt  used  a  flat  plate  specimen  to  record  the 
aerodynamic  pressure  and  heat  flux  at  the  center  of  the  plate  as  a  reference.  In  addition,  spherical 
pressure  and  thennal  domes  with  a  diameter  of  35.6  cm  and  the  three  HID  ratios  shown  in  Table 
2.1  were  instrumented.  Table  2.1  also  summarizes  the  freestream  conditions  p\  and  M\,  for  each 
test.  A  schematic  of  the  test  specimen  and  the  58  instrumented  locations  is  shown  in  Figure  2.3. 
For  the  purposes  of  this  study,  the  analysis  is  limited  to  the  points  along  the  centerline  parallel  to 
the  flow.  An  investigation  by  Ostoich  et  al.  [15,16]  discovered  that  the  recorded  data  at  points  1 
and  38  may  have  been  affected  by  an  uncharacterized  gap  between  the  dome  and  plate,  thus  only 
the  middle  1 1  data  points  along  the  centerline  (points  2-39)  were  considered. 
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l 

wd 

Figure  2.3.  Instrumentation  locations  on  spherical  dome  geometry  [44] 
Table  2.1.  Experimental  conditions  from  Glass  and  Hunt  tests  [44] 


Test 

M, 

Pi,  Pa 

D,  m 

H/D 

Boundary 

Layer 

Run  30 

6.60 

645.9 

0.355 

0.028 

Turbulent 

Run  31 

6.60 

648.0 

0.355 

0.013 

Turbulent 

Run  32 

6.59 

645.9 

0.355 

0.006 

Turbulent 

During  the  experiment,  thermocouples  recorded  temperature  measurements  at  20  samples  per 
second  over  approximately  5  seconds.  Using  an  initial  dome  temperature  Tw$  as  300K,  a  time- 
dependent  temperature  profile  across  the  dome  centerline  was  generated  using  the  linear  finite- 
difference  heat  transfer  relationship  between  heat  flux  and  temperature  change  across  time, 
Q=pCpx  ATJAt.  This  same  one-dimensional  heat  transfer  relationship  was  used  by  Glass  and 
Hunt  to  estimate  the  reported  heat  flux  from  the  thennocouple  readings,  where  r  is  the  dome 
thickness  of  0.00157m,  p  and  Cp  are  density  and  specific  heat  of  aluminum  (7000  series),  and  At 
is  equivalent  to  20  samples  per  second. 

The  linear  quasi-static  approximation  for  temperature  history  derivation  is  presented  in  Eq. 
(2.8),  also  demonstrating  the  presence  of  random  normal  measurement  error  with  zero  mean  and 
standard  deviation  of  0.5K.  Furthermore,  from  recognizing  a  linear  model  may  misrepresent  the 
true  underlying  physics  of  the  coupled  aerothennoelastic  problem,  an  exponential  discrepancy  is 
manufactured  to  approximately  simulate  reaching  equilibrium  temperature  («o=0.75K,  ai=0.5s"'), 
also  shown  in  Eq.  (2.8). 

QAt 

T’w.i+i  =  TWii  +  — —  +  a0  exp  (an)  +  er(0,  <xn)  (2.8) 

pCpT 
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Figure  2.4.  Estimated  transient  temperature  distributions  from  Runs  30,  31,  and  32 


The  temperature  profiles  over  5  seconds  in  Runs  30,  31,  and  32  are  shown  in  Figure  2.4. 
There  is  an  absence  of  two  Run  3 1  transient  temperature  histories  along  the  dome  due  to  missing 
heat  flux  measurements  those  points.  The  transient  aerothennal  data  is  used  in  this  investigation 
for  Bayesian  model  calibration  and  prediction  confidence  assessment  for  coupled,  time- 
dependent  aerothennal  analysis. 
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3  ERROR  QUANTIFICATION  AND  CONFIDENCE  ASSESSMENT  OF 
AEROTHERMAL  MODEL  PREDICTIONS 

This  section  analyzes  the  prediction  error  for  the  Glass  and  Hunt  experiments  [44]  using  the 
assumptions  and  results  from  Culler  et  al.  [11].  This  work  reevaluates  some  of  the  assumptions 
and  errors  that  were  observed  in  the  previous  study. 

First,  a  description  of  the  uncertain  input  parameters  in  the  experiments  and  aerodynamic 
pressure  and  heating  calculations  is  provided  with  sensitivity  analysis.  Next,  two  sets  of 
experimental  data  are  used  to  calibrate  uncertain  model  inputs  and  errors.  The  calibrated  inputs 
are  then  used  to  update  nominal  predictions  for  the  spherical  dome  experiments.  Then,  a  third 
data  set  is  used  for  validation  with  Bayesian  hypothesis  testing-based  confidence.  Finally,  a 
model  selection  study  is  performed  using  the  confidence  metric  for  different  forms  of  piston 
theory. 

3.1  Model  Input  Uncertainty  and  Sensitivity  Analysis 

Consider  the  flat  plate  specimen,  where  oblique  shock  relations  are  used  for  aerodynamic 
pressure  pf  and  Eckert’s  reference  temperature  method  for  aerodynamic  heating  Qf" .  Note  that 

for  the  flat  plate,  we  are  interested  in  the  value  at  the  center  of  the  plate,  which  corresponds  to 
location  ‘4’  in  Figure  2.1.  The  flat  plate  experiments  consisted  of  three  tests  (Runs  30,  31,  and 
32),  which  all  correspond  to  the  same  nominal  inputs  and  turbulent  boundary-layer  with  a  sharp 
leading  edge  panel  holder.  For  these  tests,  the  freestream  pressure  pi,  and  Mach  number  Mi, 
were  given  as  shown  in  Table  2.1.  In  addition,  the  output  aerodynamic  pressure  and  heat  flux 
were  measured  at  the  center  of  the  flat  plate.  However,  three  critical  pieces  of  information  were 
not  available  in  the  Glass  and  Hunt  report  [44]:  the  freestream  temperature  7),  wall  temperature 
Tw4,  and  equivalence  ratio  Req.  Therefore,  realistic  values  had  to  be  estimated  from  other  reports 
of  similar  testing.24  The  mean  freestream  and  wall  temperatures  are  assumed  to  be  220K  and 
300K,  respectively.  The  equivalence  ratio  is  also  uncertain,  but  for  the  current  investigation  a 
constant  value  of  Req  =  0.9  is  assumed. 

To  get  a  better  understanding  of  the  uncertainty  in  the  outputs  and  their  sensitivity  to  the 
inputs,  statistical  distributions  were  assumed  for  the  inputs.  Since  p\  and  Mi  were  measured,  1% 
coefficient  of  variation  (CV)  is  used  for  measurement  variability.  However,  10%  CV  is  used  for 
T\  and  Tw4  since  they  were  not  reported  and  had  to  be  assumed.  Normal  distributions  are  used  for 
all  four  random  inputs  and  their  distribution  parameters  are  shown  in  Table  3.1. 

Table  3.1.  Uncertainty  for  inputs  to  aerodynamic  pressure  and  heat  flux  calculations 


Measured 

Inputs 

Mean 

Standard 

Deviation 

Coefficient  of 
Variation 

Pi  (Pa) 

652.5 

6.525 

1% 

Mi 

6.6 

0.066 

1% 

Uncertain 

Inputs 

Mean 

Standard 

Deviation 

Coefficient  of 
Variation 

T\  (K) 

220 

22.0 

10% 

T w4  (K) 

300 

30.0 

10% 
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Local  and  global  sensitivity  analyses  are  performed  to  investigate  the  sensitivity  of 
aerodynamic  pressure  and  heating  to  the  input  variables.  For  defining  the  sensitivity  measures, 
let  Y  =  . .  ,Xn ) ,  where  Xt  is  the  measured  or  uncertain  inputs  and  Y is  the  resulting 

random  output.  The  local  sensitivity  is  calculated  as  the  difference  of  the  total  variance  var(7),  to 
the  variance  when  each  of  the  corresponding  random  variables  is  evaluated  at  their  mean  with 
the  other  inputs  remaining  random  (Eq.  (3.1))  [45].  The  greater  the  value  of  Act2  ,  the  greater  the 
importance  of  A)  on  Y.  Note  that  A7,  refers  to  being  calculated  over  all  random  variables  X, 
except  X,.  The  global  sensitivity  is  expressed  as  main  effect  sensitivity  index  A,-  and  total  effect 
sensitivity  index  .SV,  shown  in  Eqs.  (3.2)  and  (3.3),  respectively.21  The  Sj  of  a  variable  is  another 
measure  of  the  sensitivity  of  A)  on  Y  and  At/  provides  infonnation  about  the  interaction  of  A)  with 


other  variables.  The  sensitivities  for  the  initial  random  inputs  in  Table  3.1  are  shown  in  Table 

3.2. 

var  (7)  -  var^.  {y\X.  =  xt ) 

= - A - - 

var(T) 

(3.1) 

„  var^K^I*,)] 

var(T) 

(3.2) 

„  \X~S 

T'  var(T) 

(3.3) 

Table  3.2.  Local  and  global  sensitivity  for  aerodynamic  pressure  and  heat  flux  for  the  flat 
_ plate  with  initial  uncertainty _ 


Input  Variable 

sf 

c  pi 

Kf 

sf 

c  Q? 

ji 

Pi  (Pa) 

0.686 

0.684 

0.680 

0.020 

0.016 

0.016 

Mi 

0.333 

0.319 

0.319 

0.122 

0.174 

0.175 

T\  (K) 

0.0008 

0.0002 

0.0002 

0.451 

0.464 

0.465 

Tw4  (K) 

- 

- 

- 

0.340 

0.344 

0.344 

As  expected,  the  temperatures  play  a  small  role  in  the  pf  calculation;  Tn 4  does  not  appear  in 

oblique  shock  relations  and  T\  is  only  used  with  Req  to  determine  the  methane-air  properties. 
However,  7)  and  Tw 4  are  dominant  in  the  heat  flux  calculation  with  0.451  and  0.340, 
respectively.  Furthermore,  since  the  sum  of  the  main  effect  indices  A,  is  close  to  1,  individual 
values  of  the  main  effect  indices  S,  and  the  total  effect  indices  At,  are  so  similar,  it  is  indicated 
that  there  is  not  a  strong  interaction  among  variables.  Table  3.3  shows  the  forward  uncertainty 
propagation  of  the  normal  random  variables  from  Table  3.1  to  pf  and  Qf  . 
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Table  3.3.  Uncertainty  propagation  using  initial  uncertainty  to  pressure  and  heat  flux  for 

_ the  flat  plate _ 


Output 

Mean 

Coefficient  of 
Variation 

p?(Pa) 

1385.4 

1.21% 

Qf(  W/cm2) 

5.211 

11.53% 

Observe  that  the  10%  uncertainty  in  the  temperatures  (7)  and  7%)  play  a  larger  role  in  the 
Qf  calculation,  therefore  it  has  a  larger  CV  at  11.53%.  Since  these  experimental  values  are 
unknown  and  the  distributions  are  assumed,  it  is  beneficial  to  calibrate  these  uncertain  model 
inputs.  The  next  section  uses  Bayesian  updating  to  calibrate  the  7/  and  Tw 4  distributions  and 
quantify  the  model  errors  using  a  Bayesian  network  with  the  Glass  and  Hunt  data. 

3.2  Bayesian  Model  Parameter  Calibration 

There  is  significant  epistemic  uncertainty  in  the  true  values  of  T\  and  Tw 4,  therefore  Bayesian 
model  parameter  calibration  can  assist  in  better  approximating  these  values  based  on 
observations.  Furthennore,  the  errors  in  aerodynamic  pressure  and  heat  flux  for  the  flat  plate  and 
spherical  dome  predictions  can  also  be  calibrated.  First,  as  a  brief  introduction  to  Bayesian 
concepts,  let  (p  be  the  uncertain  model  parameters  or  errors  in  a  model  x  ( cp )  with  some  prior 

information  on  the  parameters’  uncertainty  as  a  basis  for  a  statistical  distribution  n{(p).  Then 

using  some  observed  data  y,  the  distribution  of  the  unknown  parameters  is  updated  using  Bayes 
theorem,  as  shown  in  Eq.  (3.4)  [18]. 


»(,i 

\v*\_y\x((pf\n(<p)d<p 


(3.4) 


Thus,  this  Bayesian  updating  reduces  the  uncertainty  in  the  parameters  (p ,  given  observations 
v.  In  this  case,  the  uncertain  parameters  are  (p  =  J ;  where, 

ePA  ’  ea  ’  K'i  ’  anc*  ea  are  ^1C  errors  in  the  aerodynamic  pressure  and  heat  flux  predictions  for  the 

flat  plate  and  spherical  dome  specimens.  The  model  error  is  defined  as  the  difference  between 
the  model  prediction  and  the  true  value,  as  shown  in  Eqs.  (3.5)  and  (3.6)  for  the  flat  plate  and 
spherical  dome  models,  respectively.  For  this  study,  a  systematic  error  in  the  predictions  across 
the  dome  is  assumed  for  convenience,  as  seen  in  Eq.  (3.6). 

<=<+<  &t=e.V4'  <3-5> 

J,+<  (efirfe tV<  <3-6> 
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Building  upon  the  relationship  of  the  inputs  and  predictions,  a  Bayesian  network  of  the 
measured  inputs  (p i  and  Mi),  uncertain  inputs  (7)  and  Tw4),  model  predictions  (p4  and  (Q4),  and 
model  errors  (ep  and  eQ  )  is  constructed.  Figure  3.1  depicts  the  Bayesian  network  for  the 

aerodynamic  pressure  and  heat  flux  predictions  for  the  flat  plate  and  spherical  dome  geometries 
and  the  interconnections  between  inputs,  errors,  and  data. 


Figure  3.1. 


Bayesian  network  for  calibrating  model  inputs  and  errors  using  aerothermal 

data 


The  gray  nodes  in  Figure  3.1  are  the  uncertain  inputs  and  errors  that  are  being  calibrated  with  the 
Glass  and  Hunt  data.  The  dashed-box  around  the  network  represents  the  randomness  in  the 
measured  inputs p\  and  Mi.  Bayes  theorem  in  Eq.  (3.7)  is  rewritten  for  the  case  corresponding  to 
the  aerodynamic  pressure  and  heat  flux  for  the  Glass  and  Hunt  experiments  in  Eq.  (3.7). 


Pr 

\y  p^y Q^y  p^y Mi 

\p4{(p),Q4{(p),Pi,Mx 

M  <p) 

i  y  ’ +&  ’  y  px  ’  y  Mi )  i 

H 

y  Pi  ifg,)  y  Pl  ■>  yMl  1 

p4{(p),Q4{(p),px,M^\ 

n{(p)d(p 

In  Eq.  (3.7),  uncertain  inputs  and  errors  are  (p  =  ^ 7] ,  Tv4 , J  ,  data  is  available  for  yP4 
and  yQ  from  the  flat  plate  and  spherical  dome  measurements,  as  well  as  measured  input  data  y 
and  yM  (Table  3.1).  Therefore,  all  four  sources  of  data  are  incorporated  in  the  likelihood 
function.  Now  that  (p  and  y  are  identified,  we  must  now  define  the  prior  distributions  n {(p) . 

Normal  distributions  are  used  for  uncertain  inputs  T\  and  Tw4,  with  means  from  Culler  et  al.  and 
10%  coefficient  of  variation,  as  summarized  in  Table  3.1.  Regarding  model  errors,  observations 
from  previous  reports  indicated  that  p4  and  Q4  predictions  are  expected  to  be  accurate  within 
[-10%, +  10%]  and  [- 10%, -30%] ,  respectively  [46].  The  error  bounds  for  Q4  are  associated  with 
Eckert’s  reference  temperature  method,  which  is  expected  to  consistently  over-predict  the  true 
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value  due  to  the  calorically  perfect  gas  assumption.  However,  after  a  preliminary  comparison  of 
predictions  to  data,  a  uniform  distribution  over  the  range  [-30%,  +30%]  of  the  prediction  was 
determined  to  be  a  more  appropriate  prior  for  all  four  error  terms  in  this  study.  Therefore,  this 
error  model  assumes  uniform  distributions  based  on  the  experimental  means  for  the  prior 
distribution  of  errors  x{(p)-  Normal  distributions  are  used  for  the  likelihood  function 

Pr[_y  |  x (#>)],  where  the  distribution  parameters  from  Table  3.2  are  assumed  for  yn  and  yM| , 

whereas  5%  measurement  uncertainty  is  assumed  for  aerodynamic  pressure  and  heat  flux 
measurements  yPi  and  yQt . 

Bayesian  updating  according  to  Eq.  (3.7)  is  performed  using  all  of  the  observed  data  from 
Glass  and  Hunt,  except  for  Run  30  for  the  spherical  dome.  Run  30  data  is  reserved  for  validation, 
which  is  discussed  in  the  following  section.  When  performing  the  Bayesian  updating,  the 
freestream  pressure  p\,  and  Mach  number  M\,  are  also  treated  stochastically  due  to  the 
measurement  uncertainty  presented  in  Table  3.2  with  1%  CV.  Equation  (3.7)  is  evaluated  at  100 
realizations  of  p\  and  M\  using  Latin  Hypercube  sampling.  For  each  of  those  samples,  a  Markov 
Chain  Monte  Carlo  (MCMC)  algorithm  called  slice  sampling  is  employed  using  104  samples  to 
calculate  the  posterior  distribution.  Figures  3.2  and  3.3  show  the  integrated  posterior 

distributions  for  the  uncertain  inputs  and  errors  <p  =  \^Tt,  Tw4 ,  eft , , e'^  ]  . 

The  mean  and  standard  deviation  of  the  posterior  distributions  for 
(p  =  '\Tx,  Tw4 ,  J  are  shown  in  Table  3.4.  Comparing  the  initial  and  updated 

distributions  of  T\  and  Twa,  it  is  seen  that  the  uncertainty  is  reduced,  however  the  mean  value  did 
not  shift.  This  is  primarily  a  result  of  the  errors  in  p 4  and  Qa,  predictions  being  more  easily  scaled 
as  defined  in  Eqs.  (3.5)  and  (3.6).  Thus,  calibrating  the  errors  did  result  in  a  shift  in  the  mean 
values,  as  seen  in  Figure  3.3.  Also,  there  is  significant  uncertainty  reduction  in  the  errors  from 
the  initial  ±30% . 


Figure  3.2.  Prior  and  posterior  distributions  for  a)  freestream  temperature  7),  and  b)  wall 

temperature  Tw 4 
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Figure  3.3.  Prior  and  posterior  distributions  for  a)  error  in  flat  plate  p4,  b)  error  in  flat 
plate  Q4,  c)  error  in  spherical  dome  p4,  and  d)  error  in  spherical  dome  Q4 


Table  3.4.  Mean,  standard  deviation,  and  coefficient  of  variation  of  calibrated  model 
_ inputs  and  errors _ 


Output 

Mean 

Standard 

Deviation 

Coefficient 
of  Variation 

T\  (K) 

220.67 

9.25 

4.19% 

Tw4  (K) 

300.34 

28.24 

9.40% 

<  (Pa) 

-108.87  (-8.5%) 

40.36 

37.07% 

4:  (%»■) 

-4707.0  (-10.0%) 

280.2 

5.95% 

<  (Pa) 

-75.84  (-5.0%) 

21.28 

28.06% 

<  (v/J) 

-4682.6  (-7.5%) 

247.8 

5.29% 

The  calibrated  distributions  for  T\  and  Tw4  are  propagated  to  pf  and  Qf  in  Table  3.5.  The 
uncertainty  in  aerodynamic  pressure  is  unchanged  since  it  is  insensitive  to  freestream 
temperature.  However,  the  uncertainty  in  Qf  is  reduced  from  1 1.53%  to  6.46%. 
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Table  3.5.  Uncertainty  propagation  using  updated  uncertainty  to  pressure  and  heat  flux 

_ for  the  flat  plate _ 


Output 

Mean 

Coefficient  of 
Variation 

p  f  (Pa) 

1385.4 

1.23% 

Qa  (w/cm2) 

5.229 

6.46% 

The  next  section  investigates  the  effect  of  quantifying  the  model  errors  in  the  predictions  for 
the  spherical  dome  and  uses  the  remaining  set  of  data  (Run  30)  for  assessing  the  confidence  in 
pf  and  Qf  predictions. 

3.3  Assessing  Prediction  Confidence  for  Model  Validation 

The  aerodynamic  pressure  and  heat  flux  along  the  spherical  dome  is  evaluated  at  the  initial 

and  updated  values  of  <p  =  ]  from  Table  3.6.  Figures  3.4  and  3.5  show  the 

experimental  data  from  Glass  and  Hunt  (Runs  30,  31,  and  32),  along  with  the  initial  and  updated 
model  predictions  pf  and  Qf  evaluated  at  the  mean  along  the  streamwise  centerline  of  the 

spherical  domes.  Tables  3.6  and  3.7  summarize  the  deterministic  errors  in  the  predictions.  As 
illustrated  in  Figures  3.4  and  3.5,  aerodynamic  pressure  and  heating  are  greatest  near  the  leading 
edge  of  the  dome  and  lowest  near  the  trailing  edge.  This  is  a  result  of  the  slope  of  the  dome  in  the 
flow  direction,  where  positive  slope  results  in  elevated  values  and  negative  slope  produces  lower 
values  relative  to  the  flat  plate.  Thus,  the  largest  dome  (Run  30)  produces  the  greatest  spatial 
variations  in  pressure  and  heating.  Note  that  the  slope  of  each  dome  is  zero  at  x/D= 0.5.  At  this 
location,  pressure  and  heating  values  are  nearly  identical  for  each  dome  and  for  the  flat  plate, 
which  indicates  that  local  surface  inclination  has  a  strong  impact  on  local  pressure  and  heating 
values. 


Figure  3.4.  Aerodynamic  pressure  along  centerline  of  spherical  dome  for  Runs  30,  31,  and 
32  from  test  data,  initial  mean  input  values,  and  Bayesian  updated  mean  input  values 
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Table  3.6.  Error  summary  for  nominal  pressure  predictions  along  centerline  of  spherical 

dome 


Initial  Errors  in  pf 

Updated  Errors  in  pf 

Run  30 

Run  31 

Run  32 

Run  30 

Run  31 

Run  32 

Average 

13.6% 

8.2% 

6.7% 

16.0% 

8.3% 

2.7% 

Maximum 

39.3% 

23.5% 

14.2% 

35.4% 

18.5% 

8.9% 

From  Figure  3.4  and  Table  3.6  it  is  evident  that  3ld-order  piston  theory  predictions  of  pf 
become  less  accurate  with  increasing  dome  surface  inclination.  Accordingly,  the  largest  error  in 
pf  occurs  at  the  forward-most  location  in  Run  30.  Recall  that  Run  30  was  saved  for  validation 
and  only  Runs  3 1  and  32  were  included  in  calibration.  This  generally  resulted  in  smaller  errors 
for  Runs  31  and  32,  but  errors  for  Run  30  increased,  as  summarized  in  Table  3.6.  It  is  expected 
that  if  data  from  Run  30  had  been  included  in  calibration,  then  the  corresponding  errors  would 
have  also  been  reduced.  Furthermore,  since  the  errors  in  pf  along  the  dome  vary  in  magnitude 

spatially,  it  would  be  beneficial  to  use  a  more  flexible  error  model,  such  as  a  Gaussian  process 
model,  in  more  practical  applications. 

Figure  3.5  and  Table  3.7  show  that  Eckert’s  reference  temperature  method  predicted  Qf 
more  accurately  than  3ld-order  piston  theory  predicted  pf .  Again,  the  trend  is  observed  that 
errors  are  reduced  using  the  updated  (p  values  for  Runs  31  and  32,  but  not  Run  30.  This  may  be 
expected  since  values  of  ys£  from  Runs  31  and  32  are  included  in  the  Bayesian  calibration, 
whereas  values  from  Run  30  were  not  used. 


Qf 

(w/cm1)  M 


0.0 


0.00  0.10  0.20  0.30  0.40  0.50  0.60  0.70  0.80  0.90  1.00 

x/D 

Figure  3.5.  Aerodynamic  heat  flux  along  centerline  of  spherical  dome  for  Runs  30,  31,  and 
32  from  experimental  data,  initial  mean  input  values,  and  Bayesian  updated  mean  input 

values 
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Table  3.7.  Errors  for  nominal  heat  flux  predictions  along  centerline  of  spherical  dome 


Initial  Errors  -  Qf 

Updated  Errors  -  Qf 

Run  30 

Run  31 

Run  32 

Run  30 

Run  31 

Run  32 

Average 

3.9% 

6.2% 

12.8% 

11.7% 

4.1% 

3.6% 

Maximum 

12.9% 

10.4% 

21.2% 

24.7% 

8.7% 

12.4% 

The  deterministic  errors  are  useful  for  assessing  the  accuracy  of  the  nominal  model 
predictions,  however  this  is  a  stochastic  problem  and  error  alone  does  not  provide  a  statistical 
assessment  of  the  confidence  in  the  model  prediction.  Therefore,  the  most  important  step  in  this 
model  uncertainty  framework  is  to  validate  the  models  by  assessing  the  confidence.  This  enables 
decision-making  in  regard  to  model  development  and  fidelity  selection.  For  the  Aircraft  Digital 
Twin,  it  is  important  to  have  this  confidence  metric  to  make  autonomous  decision  making 
possible  for  efficient  simulations  and  risk  mitigation. 

Several  validation  metrics  exist  with  advantages  and  disadvantages,  such  as  classical 
hypothesis  testing,  and  difference  and  area  metrics;  however  Bayesian  hypothesis  testing  is 
selected  for  this  study  [47,48],  The  Bayes  factor  approach  fits  appropriately  with  the  Bayesian 
network  integration  framework,  but  its  main  advantages  are  that  it  takes  into  account  the  entire 
probability  distribution  of  the  model  output  and  its  relation  to  a  confidence  metric  is 
straightforward.  For  Bayesian  hypothesis  testing,  we  want  to  determine  the  probability  of  our 
model  being  correct,  given  some  observed  data.  Consider  a  hypothesis  test  to  determine  the 
probability  that  a  model  prediction  x  is  equal  to  its  true  value  xq.  Equation  (3.8)  calculates  the 
Bayes  factor  B ,  as  the  ratio  of  likelihoods  corresponding  to  the  null  hypothesis  (model  prediction 
is  equal  to  the  true  value)  and  the  alternate  hypothesis  (model  prediction  is  not  equal  to  the  true 
value).  Therefore,  when  B  >  1,  the  data  supports  the  null  hypothesis  better  than  the  alternative 
hypothesis.  The  integral  fonn  of  the  Bayes  factor  in  Eq.  (3.8)  includes  the  likelihood  function  of 
the  data  supporting  the  prediction  Pr(y|x),  the  probability  density  function  (PDF)  of  the  model 
prediction  7r0  (x) ,  and  the  PDF  for  the  alternative  hypothesis  nx  (x) . 

Prfvl  H  -x  =  x)  [ Pr(  v  I  xWn  (x]dx 
B(x0)  =  y1  0  °M  u  '  ’  oK  ’  (3.8) 

V*(y\Hx\x^x0)  JPr(v|x)^1(x)Jx 

Equation  (3.8)  is  rewritten  in  Eqs.  (3.9)  and  (3.10)  for  the  cases  for  aerodynamic  pressure  and 
heat  flux,  where  /  =  1  to  1 1  for  the  x-location  across  the  dome. 
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The  Bayes  factor  is  calculated  at  each  of  the  1 1  points  along  the  spherical  dome  used  in  Run  30 
of  the  Glass  and  Hunt  study.  The  likelihood  function  Pr(y|x)  is  based  on  the  assumption  of  a 
nonnal  distribution  for  measurement  error,  with  a  standard  deviation  based  on  5%  coefficient  of 
variation  on  the  mean  of  the  pressure  and  heat  flux  data  for  Run  30  ( A(y;,0.05y) ).  The 

probability  density  function  for  the  null  hypothesis  is  detennined  by  propagating  the 

uncertainty  in  pi,  M\,  T\,  and  Tw4,  as  well  as  the  quantified  errors  for  the  spherical  dome  ( ef  and 
ef  ).  The  PDF  for  the  alternative  hypothesis  nx  (x) ,  is  modeled  as  a  unifonn  distribution 
extending  beyond  the  maximum  and  minimum  values  of  pf  and  Qsf  predictions. 

The  Bayes  factors  computed  in  Eqs.  (3.9)  and  (3.10)  can  be  used  to  the  confidence  C,  in  the 

25 

prediction,  as  shown  in  Eq.  (3.11).' 


c=-j^i=p(H«\y)  <311> 

As  indicated  in  Eq.  (3.11),  C  is  simply  the  posterior  probability  of  the  null  hypothesis  being  true, 
given  the  observation  data  (under  the  assumption  that  prior  probabilities  of  the  null  and 
alternative  hypotheses  are  both  0.5).  For  a  Bayes  factor  of  1.0,  the  confidence  C,  is  equal  to  50%. 
This  implies  that  we  do  not  have  enough  evidence  to  reject  or  accept  the  null  hypothesis. 
However,  for  Bayes  factors  greater  than  1.0  (as  explained  for  Eq.  (3.8)),  we  would  have 
increasing  confidence  that  the  prediction  is  equal  to  the  true  value.  The  confidence  metric  can  be 
used  as  a  resource  allocation  measure  for  determining  when  it  is  beneficial  to  perfonn  tests, 
where  higher  fidelity  models  are  required,  and  which  disciplines  need  a  more  strongly  (or  less 
strongly)  coupled  solution  procedure.  In  addition,  the  Bayes  factor-based  confidence  can  be  used 
to  assess  the  limits  of  the  model’s  predictions.  Table  3.8  summarizes  the  confidence  (Eqs.  (3.9)- 
(3.11))  and  detenninistic  error  in  pf  and  Qf  predictions  for  the  1 1  observations  from  Run  30  ( 
ypt  and  ). 

Table  3.8.  Error  and  confidence  in  aerodynamic  pressure  and  heat  flux  predictions  along 


centerline  of  spherical  dome  i 

or  Run  3C 

Location 

x/D 

% error 

Pt 

r 

p: t 

% error 

Qf 

r 

Ql 

1 

0.11 

-35.4% 

0.00% 

2.5% 

86.7% 

2 

0.19 

-21.0% 

0.02% 

2.8% 

87.5% 

3 

0.26 

-9.3% 

71.1% 

3.9% 

87.7% 

4 

0.34 

-3.9% 

94.0% 

4.4% 

88.6% 

5 

0.42 

2.4% 

94.6% 

8.4% 

85.1% 

6 

0.50 

7.0% 

88.7% 

9.7% 

85.0% 

7 

0.58 

11.3% 

74.8% 

14.5% 

74.1% 

8 

0.66 

15.3% 

54.9% 

21.3% 

36.1% 

9 

0.74 

21.4% 

16.4% 

24.7% 

25.0% 

10 

0.81 

23.7% 

21.7% 

17.1% 

81.8% 

11 

0.89 

25.9% 

30.2% 

19.7% 

81.5% 
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The  majority  of  the  predictions  have  greater  than  50%  confidence,  which  means  the  data 
supports  the  prediction.  Pressure  predictions  at  locations  4-6  have  the  highest  confidence.  The 
confidence  in  Qsf  predictions  using  Eckert’s  reference  temperature  method  are  all  above  80%, 

with  the  exception  of  locations  7-9.  As  mentioned  for  Figure  3.4,  the  largest  error  in  the  pf 

occurs  at  the  front  of  the  dome,  which  corresponds  to  0%  confidence.  The  deterministic  errors 
give  an  indication  of  the  quality  of  the  nominal  predictions,  however  note  that  it  is  not  always 
indicative  of  the  statistical  confidence  in  the  predicted  values.  For  example,  2.4%  error  at 
location  5  for  pf  corresponds  to  94.6%  confidence,  but  2.5%  error  at  location  1  for  Qf  is 
lower  at  86.7%  confidence.  This  is  the  result  of  differences  in  the  shape  of  the  model  error 
distributions  for  pf  and  Qf  . 

The  scope  of  this  work  is  to  initiate  a  framework  for  assessing  the  confidence  in  coupled 
aerothennoelastic  model  predictions,  not  to  necessarily  draw  definitive  conclusions  for  these 
particular  aerodynamic  models.  However,  given  the  confidence  associated  with  the  pf 

predictions  for  Run  30  {HID  =  0.028),  one  would  likely  conclude  that  3rd-order  piston  theory  is 
inadequate  for  predicting  the  aerodynamic  pressure  on  a  spherical  dome  protuberance  of  this 
size.  Although,  to  truly  reach  that  conclusion,  more  thorough  uncertainty  quantification  is 
required  to  better  capture  the  model  error. 

The  final  section  builds  upon  the  conclusions  using  the  Bayes  factor-based  confidence  and 
uses  that  metric  to  compare  predictions  using  different  fonns  of  piston  theory. 


3.4  Model  Selection  from  Prediction  Confidence  Metric 

The  Bayes  factor-based  confidence  metric  is  also  useful  for  model  selection.  Up  to  this  point, 
a  3ld-order  expansion  of  piston  theory  from  Eq.  (2.5)  was  used  to  predict  pf  .  Naturally,  1st-  and 

2nd-order  expansions  could  have  been  used  instead.  Consider  1st-  and  2nd-order  piston  theories  for 
a  model  selection  study,  as  shown  in  Eqs.  (3.12)  and  (3.13),  respectively. 


pf  =A  +  2-£- 
m3 


dw  dw 
U3  dt  dx 


(3.12) 


pf  =p,+  2-£- 

AC 


(  1  dw 
y  U3  dt  dx 


Y  + 1  ,  , 
+  - - AC 


(  1  dw  ^ 
yU3  dt  dx  j 


(3.13) 


Figures  3.6  and  3.7  show  the  pf  and  Qf  predictions  for  Run  30  using  1st-,  2nd-,  and  3rd- 
order  piston  theories.  Recall  that  Eckert’s  reference  temperature  method  uses  the  pf  predictions 
from  piston  theory,  so  Qf  is  also  affected. 
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Figure  3.6.  Aerodynamic  pressure  predictions  for  Run  30  using  1st-, 2nd-,  and  3lfl-order 

piston  theory 
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Figure  3.7.  Aerodynamic  heat  flux  predictions  for  Run  30  using  1st-, 2nd-,  and  3ld-order 

piston  theory 

Comparing  the  different  piston  theories,  both  2nd-  and  3ld-order  capture  some  of  the 
nonlinearity  in  pf  and  Qf  and  give  relatively  similar  predictions.  First-order  piston  theory 
appears  to  give  a  more  accurate  prediction  at  the  front  of  the  dome  for  pf ,  however  the  aft 
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portion  has  larger  errors  than  2nd-  and  3rd-order.  Furthermore,  the  lst-order  expansion  decreases 
the  overall  accuracy  over  the  entire  dome  for  heat  flux  predictions. 

Table  3.9  shows  the  Bayesian  hypothesis  testing-based  confidence  metric  for  pf  and  Qf 
using  1st-,  2nd-,  and  3rd-order  piston  theory.  As  expected  from  Figure  3.7,  lst-order  piston  theory 
has  higher  confidence  at  the  front  of  the  pressure  dome,  but  is  significantly  lower  when 
compared  to  2nd-  and  3ld-order  predictions  along  the  centerline  of  the  dome.  When  comparing  the 
confidence  in  2nd-  and  3rd-order  pf  and  Qf ,  not  only  are  the  two  predictions  very  similar,  but 
the  confidence  in  the  2nd-order  model  is  actually  higher  at  locations  8-11.  Therefore,  not  only 
would  a  lower-order  model  theoretically  require  lower  computational  costs,  but  it  is  statistically 
more  representative  of  observations  along  the  dome  for  Run  30. 


Table  3.9.  Confidence  in  aerodynamic  pressure  and  heat  flux  predictions  using  1st-,  2nd-, 


and  3rd-order  piston  theory  along  centerline  of  spherical  dome  for  Run  30 


Location 

x/D 

c 

Pm 

r 

lst-order 

2nd-order 

3ld-order 

lst-order 

2nd-order 

3rd-order 

1 

0.11 

0.36% 

0.00% 

0.00% 

31.5% 

84.8% 

86.7% 

2 

0.19 

63.1% 

0.09% 

0.02% 

62.2% 

86.7% 

87.5% 

3 

0.26 

94.4% 

75.6% 

71.1% 

77.8% 

87.5% 

87.7% 

4 

0.34 

95.4% 

94.0% 

94.0% 

85.7% 

88.5% 

88.6% 

5 

0.42 

94.2% 

94.7% 

94.6% 

83.8% 

85.2% 

85.1% 

6 

0.50 

89.0% 

89.1% 

88.7% 

84.7% 

85.0% 

85.0% 

7 

0.58 

68.6% 

74.9% 

74.8% 

69.3% 

73.8% 

74.1% 

8 

0.66 

20.3% 

57.1% 

54.9% 

10.6% 

38.3% 

36.1% 

9 

0.74 

0.16% 

23.3% 

16.4% 

0.25% 

35.0% 

25.0% 

10 

0.81 

0.00% 

42.7% 

21.7% 

0.23% 

88.8% 

81.8% 

11 

0.89 

0.00% 

69.0% 

30.2% 

0.00% 

92.0% 

81.5% 

3.5  Summary 

A  framework  to  quantify  the  model  error  and  assess  the  confidence  in  model  predictions  for  a 
coupled  aerothennoelastic  panel  is  outlined.  Bayesian  model  calibration,  error  quantification, 
and  prediction  confidence  assessment  procedures  are  described  for  aerothennal  models  with  data 
available  from  tests  perfonned  in  a  High-Temperature  Tunnel  on  spherical  dome  protuberances 
subjected  to  hypersonic  flow.  The  models  include  3rd-order  expansion  of  piston  theory  and 
Eckert’s  reference  temperature  method  to  predict  aerodynamic  pressure  and  heat  flux, 
respectively.  This  research  aims  to  logically  and  optimally  use  the  limited  data  available  for 
model  validation  and  decision-making.  The  freestream  and  wall  temperatures  are  assumed  since 
their  values  were  not  reported  in  the  experiments.  Bayesian  calibration  is  employed  to  update  the 
uncertain  inputs  and  quantify  errors  associated  with  aerodynamic  pressure  and  heat  flux 
predictions.  The  calibrated  input  distributions  and  quantified  model  errors  are  used  to  update  the 
model  predictions  along  the  centerline  of  a  spherical  dome  specimen.  The  infonnation  on  the 
model  error  is  used  to  calculate  the  Bayesian  hypothesis  testing-based  confidence  to  enable 
model  validation  and  model  selection  for  this  aerothennal  problem.  For  this  model  selection 
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study  among  piston  theories,  it  was  observed  that  the  highest-order  model  (3rd-order)  did  not 
result  in  the  highest  prediction  confidence  metric. 
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4  SEGMENTED  BAYESIAN  MODEL  CALIBRATION  OF  AEROTHERMAL 
MODELS 

Bayesian  updating  procedures  require  many  evaluations  of  the  models  in  the  Bayesian 
network,  which  can  be  computationally  intensive.  Most  Bayesian  calibration  work  has  been  done 
with  a  single  model  or  multiple  models  sharing  common  calibration  parameters.  Traditionally,  as 
additional  experimental  data  becomes  available  corresponding  to  a  particular  model  prediction, 
recalibration  is  perfonned  over  the  entire  network.  However,  for  a  complex  system  of  models, 
such  as  coupled  aerothermoelastic  response,  it  becomes  imperative  to  assess  the  necessity  of 
recalibrating  these  models  simultaneously  over  the  entire  network  in  contrast  to  updating  only 
the  nodes  affected  by  the  new  data  in  an  isolated,  or  segmented,  manner.  In  the  case  of  the 
Bayesian  network  in  this  research,  it  is  possible  to  use  the  data  at  the  individual  nodes  to  update 
only  the  statistically  relevant  parameters  at  that  node.  Thus,  a  segmented  Bayesian  model 
calibration  is  investigated  as  a  simplified  alternative  to  calibrating  all  of  the  models  in  the 
complex  system  simultaneously. 

4,1  Bayesian  Model  Calibration  Methodology 

In  this  section,  the  framework  to  calibrate  uncertain  model  inputs,  model  discrepancy,  and 
measurement  errors  based  on  the  observed  values  of  aerodynamic  pressure  and  heating  is 
presented. 

Consider  Bayesian  calibration  of  a  general  model  >]((/)),  which  contains  an  uncertain 
parameters  </),  with  assumed  prior  distributions  n ((/>),  based  on  prior  experience,  expert  opinion, 
or  historical  data.  With  the  observed  data  y,  the  distribution  of  the  unknown  parameters  is 
updated  using  Bayes'  theorem  in  Eq.  (4.1). 

*(t\  (4.1) 

|  Pr[y  I 

Thus,  Bayesian  updating  reduces  the  uncertainty  in  parameters  (j)  given  observations  y.  Bayes’ 
theorem  is  rewritten  for  the  aerodynamic  pressure  and  heat  flux  predictions  and  Glass  and  Hunt 
experiments  in  Eq.  (4.2). 

,  , ,  ,  Pfiyp  >  To  I  Po  s  (A  Pvt  (A  Qert 

I  Tp>Tq)  =  7 - ± -  (4.2) 

J  pr[yp,yQ  I  Posit), Prr(f)>QiiKr(t)Mf)dj 

The  Bayesian  network  for  Eq.  (4.2)  is  constructed  using  the  model  and  data  relationships  of 
the  aerothennal  models  presented  in  Section  2.1.  As  previously  mentioned,  oblique  shock 
relations  predict  p 3  at  the  leading  edge  of  the  test  specimen  (i.e.  flat  plate  or  dome  protuberance). 
For  the  flat  plate,  the  pressure  across  the  panel  is  the  same  as  the  leading  edge  pressure,  therefore 
p/p  =  P  fp  ■  Thus,  the  flat  plate  pressure  data  can  be  directly  applied  to  the  oblique  shock 
prediction  such  that  yf  =  /?os(<A)  +  £os  +  N(0,av  ) .  This  isolates  the  known  relationship  of  the 

p  yp 

oblique  shock  prediction  to  the  flat  plate,  reserving  the  rigid  dome  data  to  calibrate  other 
components  of  the  aerothennal  model.  Similarly,  piston  theory  and  Eckert’s  reference 
temperature  models  can  be  related  to  Glass  and  Hunt  data.  In  effect,  pressure  data  across  the 
spherical  dome  ( y s* )  is  applied  to  calibrate  the  model  discrepancy  associated  with  piston  theory 

(£pt)  and  heat  flux  data  (y  p ,  ySg  )  is  used  to  calibrate  the  discrepancy  associated  with  Eckert’s 

reference  temperature  method  (£ERT  )•  Note  that  this  description  of  separating  the  models  and 
data  is  used  to  motivate  and  explain  the  segmented  calibration  approach,  however  simultaneous 
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calibration  combines  all  of  the  model  predictions  and  data  into  a  single  calculation.  The 
segmented  model  calibration  procedure  is  further  discussed  in  the  next  section. 

Figure  4.1  shows  the  Bayesian  network  and  the  relationships  between  the  aerodynamic 
pressure  and  heat  flux  model  predictions  (p.t,  p4,  Q4),  Glass  and  Hunt  data  ( yf ,  ys‘‘ ,  }p  ,  y'‘‘ ), 
model  inputs  (pi ,  Mi)  and  uncertain  model  inputs  (T t,  Tw  ),  measurement  errors  ( sy  ,£yg),  and 
discrepancy  terms  for  calibration  (sos,  £pt,  £ert). 


Figure  4.1.  Bayesian  Network  for  Glass  and  Hunt  experiments 

Recall  p4  prediction  from  piston  theory  is  a  function  of  the  oblique  shock  results  p 3  in  Eq.  (2.5). 
It  is  then  essential  to  consider  the  propagation  of  the  error  in  oblique  shock  80s  in  the  P4 
calculation  as  shown  in  Figure  4.1.  Similarly,  Eckert's  reference  temperature  uses  the  T4  obtained 
from  isentropic  relations  with p4,  warranting  the  passing  of  Sp-pto  the  heat  flux  model. 

4.2  Model  Discrepancy  Formulation 

The  calibration  framework  followed  in  this  study  is  based  on  the  Kennedy  and  O’Hagan 
discrepancy  function  [21],  This  framework  has  the  capacity  to  account  for  the  relevant  sources  of 
uncertainty  identified  in  this  study  and  include  them  in  the  calibration.  Their  combined  form  is: 

y  =  p(x,</>)  +  S(x)  +  £  (4.3) 

where,  ;/  is  the  model  under  consideration,  (/)  are  the  parameters  being  calibrated,  S  is  the  model 
inadequacy  function  evaluated  at  a  specific  inputs  x  and  s  \  A,r(0,s)is  the  measurement 
uncertainty  associated  with  that  set  of  data.  In  line  with  the  discussion  of  relevant  sources  of 
uncertainty  (i.e.  data  uncertainty,  input  parameter  uncertainty,  input-specific  model  inadequacy) 
the  parameters  calibrated  using  the  Glass  and  Hunt  hypersonic  wind  tunnel  data  are  listed  in  Eq. 
(4.4). 


^  5  Gts  ’  G>t  ’  Grt  ’  G 


yQ 


(4.4) 


where  eos,  £pt,  and  £ert  are  discrepancy  models  for  oblique  shock  relations,  piston  theory,  and 
Eckert’s  reference  temperature  method,  respectively.  The  measurement  errors  sy  and  sY  reflect 
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the  uncertainty  in  aerodynamic  pressure  and  heating  observations  (yp  and  yo)  observations, 
respectively. 

In  a  related  study  by  Smarslok  et  al.  [23],  systematic  error  models  were  chosen  to  represent 
the  inadequacy  in  the  aerodynamic  pressure  and  heat  flux  predictions.  It  was  observed  that  these 
error  models  were  not  fully  sufficient  to  predict  the  Glass  and  Hunt  data  with  reasonable 
confidence  across  the  dome.  Thus,  a  higher  order  model  for  discrepancy  is  investigated  for  piston 
theory  and  Eckert’s  reference  temperature  predictions. 

Figures  4.2  and  4.3  show  the  nominal  aerothermal  model  predictions  compared  to  Glass  and 
Hunt  data  for  Runs  30,  31,  and  32.  The  nominal  values  of  T\  and  Tw  are  220  K  and  300  K, 
respectively. 


3500 


3000 


-- 

-Run  30  Pred. 

Run  31  Pred. 

— 

-  Run  32  Pred. 

• 

Run  30  Data 

■ 

Run  31  Data 

♦ 

Run  32  Data 

Figure  4.2.  Aerodynamic  pressure  predictions  across  dome  at  nominal  input  values 


Figure  4.3.  Aerodynamic  heat  flux  predictions  across  dome  at  nominal  input  values 
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Table  4.1.  Errors  at  nominal  aerothermal  model  predictions  and  Glass  and  Hunt  data 


Prediction 

Run 

Maximum 

Error 

Maximum 
%  Error 

Aerodynamic 

Run 

1149 

59.3% 

Pressure 

30 

Run 

522 

34.1% 

31 

Run 

253 

17.8% 

32 

Aerodynamic  Heat 

Run 

-6549 

-14.7% 

Flux 

30 

Run 

9557 

15.2% 

31 

Run 

8210 

15.3% 

32 

Recalling  that  piston  theory  (Eq.(2.5))  is  a  function  of  slope  (fiw/fix)  and  the  pressure 
prediction  error  is  proportional  to  the  dome  slope,  we  consider  discrepancy  models  for  piston 
theory  and  Eckert's  reference  temperature  that  are  also  functions  of  slope.  In  adherence  to  the 
Kennedy  and  O’Hagan  framework,  the  uncertain  parameters  in  Eq.  (4.4)  have  been  revised  to 
reflect  the  slope  dependency  of  epj  and  £ert  in  Eqs.  (4.5)  and  (4.6),  respectively. 


=  60PT  +  6PT— +  ZyPT(— )2 
OX  ox 


ERT  .  ERT 
^ERT  ~  CQ  C\ 


dw  ert.Sw  2 

—  +  c,  ( — ) 
fix  fix 


(4.5) 

(4.6) 


Quadratic  models  are  chosen  to  represent  the  relationship  between  spherical  dome  slope  and 
Pa  and  Qa  prediction  discrepancies,  where  the  quadratic  coefficients  are  calibrated.  Also,  defining 
measurement  variability  as  sy  ~  N( 0,av  )and  s  ~  N(0,cr  )  ,  the  final  form  of  the  calibration 

parameters  is  presented  in  Eq.  (4.7). 


<t>  = 


T  T  £  £ 

Jl’%’eOS’ePT 


'  dw 

|bPT^ 

f  3W  CERT  'l 

,  G  (  CT  ) 

\,£  (cr  | 

v  fix 

J 

?  ^ERT 

Idx  ) 

yp\  yp  ) 

yQ  \  yQ  ) 

(4.7) 


4.3  Uncertain  Inputs,  Discrepancy  Parameters,  and  Measurement  Variability 

In  this  study,  normal  prior  distributions  are  used  for  uncertain  inputs  Tx  and  Tw  with  mean 

values  from  Culler  et  al.  [9]  and  10%  coefficient  of  variation  assumed  by  Smarslok  et  al.  [23],  as 
shown  in  Table  4.2.  The  measurement  uncertainty  terms  are  assumed  to  have  lognormal  prior 
distributions  (due  to  the  requirement  they  always  be  positive),  with  log-means  of  around  5%  of 
the  mean  values  of  y^andy^ .  Regarding  model  discrepancy  tenns,  uniform  distributions  are 
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assumed  with  bounds  shown  in.  These  bounds  were  obtained  by  propagating  model  input  and 
data  uncertainty  through  the  network  to  detennine  sufficient  ranges  for  the  parameters. 

Table  4.2.  Prior  Distributions  for  Uncertain  Input  Parameters 


Parameter 

Normal  jU 

cv 

Ti{  K) 

220 

10% 

Tw  (K) 

300 

10% 

Lognormal 

Lognormal 

M 

a 

% 

4.20 

0.17 

ayQ 

7.81 

0.24 

Lower 

Upper 

Bound 

Bound 

^os(Pa) 

-330 

330 

h0PT(Pa) 

-200 

200 

*T(Pa) 

-10000 

0 

hf(Pa) 

-100000 

30000 

cERT(W/m2) 

-20000 

20000 

CjERT(W/m2) 

0 

300000 

cERT(W/m2) 

2000000 

5000000 

4.4  Problem  Evaluation  for  Implementation  of  Segmented  and  Simultaneous  Calibration 
Strategies 

For  the  aerothennal  Bayesian  network  in  Figure  4.1,  the  subdivision  of  Glass  and  Hunt  data 
and  its  application  to  specific  aerothennal  models  is  the  first  critical  step  to  warrant  a  similarly 
segmented  calibration  procedure.  Intuition  and  practicality  of  segmented  calibration,  along  with 
projected  computational  savings,  prompts  the  investigation  of  its  viability  in  this  section. 

In  order  to  implement  the  segmented  calibration  strategy,  it  is  important  to  discern  which 
uncertain  variables  make  significant  contributions  to  the  overall  uncertainty  in  the  each  model 
output.  To  accomplish  this,  a  global  stochastic  sensitivity  analysis  is  performed.  Global 
sensitivity  analysis  is  based  on  the  variance  decomposition  theorem  stated  below  in  Eq.  (4.8). 
This  simply  states  that  the  variance  in  a  model  output  Y  can  be  decomposed  into:  1)  the  variance 
of  the  expected  values  of  Y  conditioned  on  a  fixed  input  quantity  (X1)  and  all  other  calibration 
quantities  (X)  allowed  to  vary;  and  2)  the  expected  value  of  the  variance  of  Y  conditioned  on  the 
same  set  [45,49]. 


V(Y)  =  V(E(Y  |  X1))  +  E{V{Y  \  X'))  (4.8) 

In  Eq.  (4.8),  V(Y)  is  computed  by  evaluating  the  total  variance  of  the  model  output 
considering  uncertainty  in  all  parameters.  The  total  effects  of  a  parameters  uncertainty  contain 
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the  interactions  between  the  7th  parameter  and  the  other  parameters.  The  interaction  between 
parameters  further  induces  uncertainty  in  the  output,  thus  a  parameter's  total  effects  are  essential 
to  gain  insight  into  the  total  significance  of  its  uncertainty.  The  total  effects  are  computed  from 
Eq.  (4.9).  For  piston  theory  and  heat  flux  prediction  sensitivities  to  the  uncertain  parameters 
across  the  dome,  an  average  total  sensitivity  value  is  reported  in  Tables  4. 3 -4. 5. 


Ex-i(Vx,(Y  \X~‘)) 
V(Y) 


(4.9) 


The  total  effect  sensitivities  are  computed  in  a  double  nested  loop  and  are  computationally 
demanding  when  calculating  heat  flux  prompting  the  construction  of  a  linear  regression  surrogate 
model  of  Eckert's  reference  temperature.  Although  some  accuracy  in  the  sensitivity  calculation 
was  lost,  the  effects  were  minimal  and  the  computational  cost  was  greatly  reduced.  The  surrogate 
heat  flux  model  used  for  sensitivity  analysis  can  be  found  in  Eq.  (4.10),  where  Qf  =  [ 1 ,  p„  Th  Mb 

7V]  and  i  is  the  location  along  the  dome.  Coefficients  corresponding  to  each  heat  flux  prediction 
point  along  the  dome  b;  were  found  using  least-squares  regression  with  200  training  points.  The 
average  R ~  value  over  the  dome  for  the  Eckert’s  reference  temperature  surrogate  model  was 
0.998. 


Qi  =  Qfb, 


(4.10) 


Tables  4. 3-4. 5  show  the  total  effect  sensitivities  using  the  prior  distributions  for  each  of  the 
uncertain  parameters  and  errors  for  the  models. 

Table  4.3.  Oblique  shock  (p3 )  sensitivities  to  prior  distributions 


Parameter  SjOS 

T,  (K)  1.5e-3 

gos(Pa)  0.999 


Table  4.4.  Piston  theory  (/j4)  sensitivities  to  prior  distributions 


Parameter  5"rpx 

T,  (K)  l.5c-3 

gos  (Pa)  0.685 

gPT  (Pa) _ 0.310 
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Table  4.5.  Eckert's  reference  temperature  ( Q4 )  sensitivities  to  prior  distributions 


Parameter 

C* 

*JT,ERT 

Ti(  K) 

0.075 

Tw  (K) 

0.007 

£os  (Pa) 

0.547 

err  (Pa) 

0.045 

eERT(W/m-) 

0.333 

From  these  sensitivities,  we  can  gauge  the  effectiveness  of  calibrating  a  particular  parameter 
with  a  particular  segmented  model  and  data  set.  For  example,  in  Table  4.3,  the  oblique  shock  and 
piston  theory  models  are  insensitive  to  T] .  Therefore,  we  would  not  expect  observations  of  yf 

and  ySp  to  infonn  us  about  the  true  value  of  T\.  It  is  also  expected  that  due  to  the  low  sensitivity 

of  all  models  to  the  prior  distributions  of  T )  and  Tw  compared  to  the  model  errors,  the  posterior 
distributions  of  uncertain  input  parameters  are  not  expected  to  update  as  much  as  those  for  the 
discrepancy  parameters.  The  segmented  calibration  procedure  will  be  conducted  in  recognition 
of  these  sensitivities  as  follows:  cr  will  be  paired  with  the  piston  theory  for  calibration,  and  T), 

y p 

Tw,  and  cr  with  Eckert’s  reference  temperature  model. 


4.5  Simultaneous  and  Segmented  Calibration  Results 

Simultaneous  and  segmented  Bayesian  model  calibrations  were  conducted  as  detailed  in  the 
previous  section  and  indicated  in  Figure  4.3.  For  simultaneous  calibration,  150,000  samples  were 
generated  using  slice  sampling18.  Similarly,  for  segmented  calibration  100,000  samples  were 
generated  for  each  segment.  The  posterior  distributions  of  the  calibration  parameters  (model 
inputs,  discrepancy  parameters,  and  measurement  errors)  are  shown  in  Figures  4. 4-4. 6. 


Prior 

- Simultaneous 

- Sequential 


Figure  4.4.  Prior  and  posterior  distributions  for  a)  freestream  temperature  and  b)  wall 

temperature 
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ER7 


0.01 


Figure  4.5.  Prior  and  posterior  distributions  for  (a)  oblique  shock  error,  (b)-(d)  piston 
theory  error  and,  (e)-(g)  Eckert's  reference  temperature  error 


Figure  4.6.  Prior  and  posterior  distributions  for  (a)  pressure  and  (b)  heat  flux 

measurement  error 
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In  Figures  4. 4-4. 6,  simultaneous  and  segmented  calibrations  produce  similar  marginal  posterior 
distributions;  the  differences  will  be  quantified  in  the  next  section  with  Bayesian  hypothesis 
testing.  The  posteriors  of  the  uncertain  inputs  (7)  and  Tw)  do  not  differ  much  from  the  prior 
distributions,  which  was  expected  due  to  the  low  sensitivities  to  those  parameters  addressed  in 
Table  4.3.  However,  the  posteriors  of  the  discrepancy  model  parameters  and  measurement  errors 
updated  significantly  from  the  priors  with  both  simultaneous  and  segmented  calibration. 

Preceding  the  calibration,  the  sensitivities  to  the  prior  parameter  distributions  were  presented 
in  Table  4.2.  The  posterior  sensitivities  are  shown  in  Tables  4. 6-4. 8  where  it  can  be  seen  that  the 
sensitivity  to  eos  has  decreased  in  piston  theory  and  Eckert’s  reference  temperature  calculations 
and  is  transferred  to  the  sensitivity  in  £PT  and  £Ert,  respectively.  Also,  because  of  the  reduced 
uncertainty  in  the  discrepancy  parameters,  the  Eckert’s  reference  temperature  model  becomes 
more  sensitive  to  T\,  however  low  sensitivity  persists  for  Tw.  The  lack  of  substantial  updating  for 
Tw  is  in  agreement  with  results  from  Smarslok  et  al.  [23]. 

Table  4.6.  Oblique  shock  (p3 )  sensitivities  to  simultaneous  calibration  posterior 

distributions 


Parameter 

O* 

JT,OS 

Ti  (K) 

0.005 

sos  (Pa) 

0.998 

Table  4.7.  Piston  theory  (p 4)  sensitivities  to  simultaneous  calibration  posterior  distributions 


Parameter 

C' 

T  PT 

T,  (K) 

0.007 

£os  (Pa) 

0.397 

ePT  (Pa) 

0.606 

Table  4.8.  Eckert's  reference  temperature  ( Q4 )  sensitivities  to  simultaneous  calibration 

posterior  distributions 


Parameter 

O* 

*JT,ERT 

T,  (K) 

0.315 

Tw  (K) 

0.040 

£os  (Pa) 

0.144 

ePT  (Pa) 

0.010 

£ert  (W/m-) 

0.471 

4.6  Comparison  of  Calibration  Strategies 

To  assess  the  effectiveness  of  the  segmented  calibration  methodology  compared  to  the 
simultaneous  procedure,  it  is  necessary  to  specify  the  comparison  criteria.  As  previously  stated, 
the  purpose  of  this  study  is  to  explore  the  computational  effort  and  accuracy  of  segmented  and 
simultaneous  calibration  methods.  Computational  effort  is  gauged  by  the  rates  of  convergence 
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between  the  two  procedures,  and  accuracy  is  determined  by  the  resulting  posterior  distributions’ 
comparative  ability  to  predict  Run  30  of  the  Glass  and  Hunt  data,  which  is  calculated  using 
Bayesian  hypothesis  testing.  Recall  that  the  calibration  was  performed  using  Runs  31  and  32, 
which  have  lower  H/D  ratios  than  the  Run  30  data,  thus  making  it  an  extrapolation  case  for 
validation. 

Calibration  convergence  rates  are  computed  using  Kullback-Leibler  (K-L)  distance  [50], 
presented  in  Eq.  (4.11).  This  metric  compares  the  final  posterior  distribution  obtained  from  the 
calibration,  7To((f>),  with  the  distribution  at  the  7th  sample,  7r,fi/>).  Convergence  is  tested  at  every 
1,000  sample  points  and  is  reached  when  Dklj  <  0.01.  The  number  of  samples  at  which  the 
simultaneous  and  segmented  procedures  reached  convergence  is  presented  in  Table  4.9. 

DklA^o\\^)=  J  x0(f  )\og^-[df  >0  (4.11) 

Table  4.9.  Segmented  and  Simultaneous  Calibration  Samples  to  Convergence 


Procedure 

Model 

Samples 

Simultaneous 

All 

73,000 

Segmented 

Oblique  Shock 
Piston  Theory 
Eckert’s  Ref. 

15,000 

25,000 

55,000 

Temp. 

Table  4.9  shows  that  in  order  to  obtain  convergence  in  simultaneous  calibration,  all  three 
models  must  be  evaluated  73,000  times.  Whereas,  for  segmented  calibration  to  obtain 
convergence  according  to  K-L  distance,  15,  25,  and  55  thousand  samples  were  required  for 
oblique  shock,  piston  theory,  and  Eckert’s  reference  temperature  method,  respectively.  The 
computational  savings  by  calibrating  with  the  segmented  approach  is  from  individual  models  not 
needing  to  be  executed  as  many  times  to  converge.  However,  the  actual  savings  from  segmented 
calibration  is  model  dependent  and  not  fully  quantified  in  this  study. 

In  Bayesian  hypothesis  testing,  we  compute  the  likelihood  ratio  of  two  competing  models 
given  observed  data.  Equation  (4.12)  calculates  this  ratio  of  likelihoods  corresponding  to  the  null 
hypothesis  (segmented  calibrated  parameters  </)<)  with  joint  probability  distribution  tto  are  equal  to 
the  true  value)  and  the  alternate  hypothesis  (simultaneous  calibrated  parameters  (j)\  with  joint 
probability  distribution  tc\  are  equal  to  the  true  value).  Therefore,  when  B(x\)  =  1  the  segmented 
calibration  prediction  has  the  same  likelihood  of  occurring  as  the  simultaneous  calibration 
prediction  at  the  data  location  x,.  The  integral  form  in  Eq.  (4.12)  includes  the  likelihood  function 
of  the  data  supporting  the  prediction  Pr(y\(f>)  the  probability  density  function  (PDF)  of  the  null 
hypothesis  7To(^o)  and  alternative  hypothesis  n\((f>\  ). 


_  Prjy  |  #,,x.)  _  J P,iy  I 

J Pr(y  |  <j\,x^)nx{<j>[)d<^ 


(4.12) 
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First,  Bayesian  hypothesis  testing  was  perfonned  separately  using  flat  plate  pressure  data, 
spherical  dome  pressure  data,  and  heat  flux  data  to  gauge  how  well  segmented  calibration 
perfonned  for  each  model.  Then,  a  Bayes  factor  was  computed  using  all  models  and  data  to 
obtain  a  combined  assessment  of  segmented  calibration.  Bayes  factors  were  obtained  at  each 
dome  location,  of  which  the  average,  minimum,  and  maximum  are  reported  for  piston  theory  and 
Eckert’s  reference  temperature  predictions  across  the  dome.  The  trends  of  these  individual  Bayes 
factors  are  values  lower  than  one  on  the  outer  dome  edges  and  greater  than  one  closer  to  the 
center.  Figures  4.7-4. 8  show  that  while  the  mean  predictions  are  very  similar,  the  uncertainty  in 
the  prediction  close  to  the  dome  center  using  the  segmented  posterior  distribution  is  much  larger 
than  the  uncertainty  using  simultaneous  posteriors,  and  thus  has  higher  likelihood  of  containing 
the  data.  This  accounts  for  the  Bayes  factors  favoring  segmented  calibration  for  predictions  near 
the  center  of  the  dome.  The  results  are  reported  in  Table  4.10. 

Table  4.10.  Bayes  Factors  for  segmented  and  simultaneous  calibration 


Model 

Minimum 

Maximum 

Average 

Oblique  Shock 

- 

- 

0.96 

Piston  Theory 

0.77 

1.19 

0.92 

Eckert’s  Ref. 

0.84 

0.81 

1.08 

0.94 

1.15 

Temp. 

Combined 

1.89 

Figure  4.7.  Pressure  predictions  across  dome  from  posterior  distributions 
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Figure  4.8.  Heat  flux  predictions  across  dome  from  posterior  distributions 


Table  4.11.  Errors  in  calibrated  aerothermal  model  predictions  and  Glass  and  Hunt  data 


Prediction 

Maximum 

Maximum 

Error 

%  Error 

Aerodynamic 

Pressure 

-191 

-10.4% 

Aerodynamic  Heating 

-12493 

-13.1% 

Table  4.11  shows  the  effect  of  the  calibration  on  the  errors  in  the  pressure  and  heating 
prediction  compared  to  the  Glass  and  Hunt  data.  The  maximum  percent  error  in  the  nominal 
aerodynamic  pressure  prediction  across  the  dome  for  Run  30  is  59.3%  which  decreased  to  - 
10.4%  with  the  calibrated  models.  Similarly,  while  the  magnitude  of  the  maximum  heat  flux 
error  for  Run  30  increased,  the  maximum  percent  error  in  the  heat  flux  prediction  decreased  to  - 
13.1%  from  -14.7%  at  the  nominal  prediction. 

Another  explanation  for  the  differences  between  predictions  using  segmented  and 
simultaneous  calibration  posteriors  is  that  segmented  calibration  ignores  correlation  between 
model  error  parameters.  These  correlations  are  significant  in  the  simultaneous  posterior  samples 
and  are  listed  in  Table  4.12.  For  example,  when  i'os  is  small  the  first  term  in  epj  is  large  in  the 
simultaneous  posterior  samples.  However,  for  predictions  from  the  segmented  calibration 
posteriors,  we  select  samples  of  eos,  £pt  ,  and  £ert  independently,  resulting  in  a  larger  variance  in 
prediction.  Table  4.12  lists  the  most  significant  correlations  between  £os,  terms  in  £PTand  terms  in 
£ert- 
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Table  4.12.  Correlations  between  model  error  parameters  in  simultaneous  posterior 
_ samples _ 


Parameter  1 

Parameter  2 

Correlation 

^os(Pa) 

60PT(Pa) 

-0.91 

^os(Pa) 

(Pa) 

-0.51 

^os(Pa) 

*TC  Pa) 

-0.13 

^os(Pa) 

c0ERT(W/m2) 

-0.35 

^os(Pa) 

CjERT(W/m2) 

-0.17 

G)s(Pa) 

cERT(W/m2) 

0.11 

C(Pa) 

cERT(W/m2) 

0.26 

C(Pa) 

CjERT(W/m2) 

0.14 

Pa) 

cERT(W/m2) 

-0.03 

4.7  Conclusions 

A  segmented  Bayesian  model  calibration  approach  is  investigated  as  an  alternative  to  full, 
simultaneous  calibration  for  isolated  uncertainty  quantification  and  reduced  computational  cost. 
The  objectives  in  this  study  of  a  segmented  calibration  approach  include  identifying  the  required 
characteristics  of  calibration,  identifying  the  appropriate  uncertain  parameters  and  errors  for 
calibration  throughout  the  segmented  process,  and  assessing  the  improved  efficiency  and 
viability  of  segmented  model  calibration.  A  Bayesian  network  is  constructed  for  aerodynamic 
pressure  and  heating  model  predictions  corresponding  to  aerothennal  experiments  for  flat  plate 
and  spherical  dome  protuberances  subjected  to  hypersonic  flow  conditions.  The  aerodynamic 
pressure  and  heat  flux  data  are  used  for  simultaneous  and  segmented  calibration  of  uncertain 
model  inputs,  measurement  errors,  and  model  discrepancy  for  aerothennal  predictions.  The 
aerothennal  problem  was  segmented  into  oblique  shock  relations,  piston  theory,  and  Eckert’s 
reference  temperature  method.  To  quantify  the  viability  and  potential  benefit  of  isolating 
calibrations  of  models  in  the  network,  segmented  and  simultaneous  calibration  are  compared 
using  the  Kullback-Leibler  distance  and  Bayes  factor  metrics.  For  model  calibration  using  the 
aerothennal  data,  the  segmented  approach  yielded  greater  prediction  uncertainty  than  the 
simultaneous  approach  due  to  inherent  conelations  lost  through  the  segmentation.  However,  the 
Bayes  factor  values  comparing  the  likelihoods  of  simultaneous  and  segmented  calibration  results 
are  still  close  to  1,  therefore  neglecting  correlations  may  be  acceptable  for  this  problem.  In 
addition,  the  reduction  in  sample  size  required  for  convergence  using  the  segmented  approach 
instead  of  simultaneous  calibration  was  79.5%,  65.8%,  and  24.7%  for  the  three  aerothennal 
models,  respectively,  needed  for  convergence  using  the  simultaneous  approach. 

Further  studies  using  the  aerodynamic  pressure  and  heat  flux  data  to  calibrate  models  in  the 
aerothermoelastic  prediction  will  be  developed  to  include  the  heat  transfer  and  aeroelastic 
disciplines..  For  example,  the  Glass  and  Hunt  data  can  be  extended  to  calibrate  the  transient  heat 
transfer  model  parameters  and  time-dependent  discrepancy  in  addition  to  the  aerothennal  models 
considered  in  this  study.  Other  investigations  of  interest  are  the  effect  of  coupling  between 
aeroelastic  and  aerothennal  predictions  as  well  as  the  dependence  the  calibration  parameters  on 
model  assumptions  and  simplifications.  Once  all  models  used  in  the  aerothennoelastic  prediction 
are  calibrated,  confidence  in  the  updated  prediction  needs  to  be  quantified. 
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5  BAYESIAN  CALIBRATION  OF  AEROTHERMAL  MODELS  USING  TIME- 
DEPENDENT  DATA 

In  this  section,  the  Bayesian  calibration  approach  for  using  observed  aerodynamic  pressure 
and  temperature  measurements  to  reduce  uncertainty  in  the  aerothermal  predictions  is  presented. 
Uncertain  model  inputs,  model  discrepancy,  and  measurement  errors  are  considered  for 
calibration.  Additionally,  the  dynamic  Bayesian  network  is  constructed  for  the  coupled,  transient 
aerothermal  problem  and  two  implementations  of  time-dependent  discrepancy  models  are 
compared  and  discussed.  Prior  distributions  for  discrepancy  model  parameters  in  Eckert’s 
reference  temperature  and  linear  heat  transfer  are  presented.  A  direct  calibration  perfonnance 
comparison  between  global  and  incremental  discrepancy  application  at  the  dome  midpoint 
follows  using  a  time-dependent  Bayes’  factor  for  model  selection.  The  selected  discrepancy 
treatment  is  applied  to  the  full  dome  for  validation. 

5,1  Bayesian  Model  Calibration  Methodology  for  Time-Dependent  Data 

Consider  Bayesian  calibration  of  a  model  ;/(</>),  which  contains  an  uncertain  parameters  (/) 
composed  of  model  inputs,  discrepancy  parameters,  or  measurement  errors  with  assumed  prior 
distributions  n(<f>),  based  on  prior  experience,  expert  opinion,  or  historical  data.  With  the 
observed  data  y,  the  distribution  of  the  unknown  parameters  is  updated  using  continuous  Bayes' 
theorem  in  Eq.  (5.1). 


y  /  Pr[y|r?(0)]7r(0)d0 

Thus,  Bayesian  updating  reduces  the  uncertainty  in  parameters  (j)  given  observations  y. 
Bayes’  theorem  is  rewritten  for  the  aerodynamic  pressure,  heat  flux,  and  heat  transfer  predictions 
(from  piston  theory,  Eckert’s  reference  temperature,  and  linear  heat  transfer,  respectively)  and 
Glass  and  Hunt  data  in  Eq.  (2.8).  Recall  from  Section  2,  the  heat  flux  data  reported  by  Glass  and 
Hunt  was  used  to  extract  the  temperature  history  across  the  panel.  For  this  reason,  only  the 
derived  temperature  history  and  pressure  can  be  used  in  the  Bayesian  network. 

/  ■  n  _  Pr[yP,yrlPpr(</>)>Q£RT(</>)>Ptfr(</>)]^(0)  (5  2) 

yP’yT  f  Pr[yp,yT\pPT(4>),QERT(4>),THT(4>)]jr(4>)d4> 

The  Bayesian  network  for  Eq.  (5.2)  is  constructed  using  the  relationships  between  the 
aerothennal  models  and  available  Glass  and  Hunt  data  presented  in  Section  2.2.  Figure  5.1 
shows  the  Bayesian  network  and  the  relationships  between  the  aerodynamic  pressure,  heat  flux, 
and  heat  transfer  models  (p4,  (A,  Tw),  aerothermal  data  (yp>  yT),  model  inputs  (pi,  M\_  TW: o), 
measurement  errors  (syT,eyp),  and  discrepancy  terms  for  calibration  (£pT,  £-ErT,  £ht)- 
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Figure  5.1.  Aerothermal  Bayesian  Network  for  Glass  and  Hunt  experiments 

Recall  that  the  Qa  prediction  from  Eckert's  reference  temperature  is  a  function  of  T\  obtained 
from  isentropic  relations  with  p\  predicted  by  piston  theory.  Therefore,  the  error  in  piston  theory 
£pi  is  propagated  to  Qa  calculation,  as  shown  in  Figure  5.1.  Similarly,  the  error  in  the  heat  flux 
boundary  condition  prediction  should  be  passed  it  forward  to  the  heat  transfer  model. 

Note,  the  Bayesian  network  in  Figure  5.1  has  both  stationary  and  transient  components,  with 
pressure  remaining  constant  for  a  rigid  dome  configuration,  while  temperature  and  heat  flux 
evolve  through  time.  The  transient  module  in  Figure  5.1  can  be  expanded  to  a  transient  Bayesian 
network  of  time-dependent  heat  flux  and  temperature  predictions  and  data  as  represented  in 
Figure  5.2.  The  transient  Bayesian  network  illustrates  the  relationship  between  the  available 
models,  Glass  and  Hunt  data,  and  model  discrepancy  e  in  a  coupled,  quasi-static  prediction. 


Figure  5.2.  Bayesian  Network  with  Time-Dependent  Data  for  Aerothermal  Coupling 

The  primary  objective  of  this  study  is  to  observe  the  coupled,  time-dependent  aerothermal 
models  and  investigate  strategies  for  modeling  the  discrepancy.  For  example,  using  these 
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models  in  succession  for  a  transient  aerothermal  analysis  without  correction  results  in  the 
coupled  aerothermal  prediction  and  temperature  data  comparison  shown  in  Figure  5.3.  Observe 
that  the  uncorrected  models  approximately  predict  a  linear  temperature  history  for  the  dome 
midpoint  and  cannot  capture  the  system’s  induced  movement  toward  equilibrium.  Figure  5.4 
demonstrates  the  discrepancy  magnitude  growing  through  time  with  a  maximum  error  of  18 
degrees  Kelvin  observed  at  5  seconds. 


Aerotherma  I  pred  ictio  n  *  M  idpo  i  nt  data 

Figure  5.3.  Nominal  aerothermal  temperature  prediction  compared  to  midpoint  data 


Time(s) 

Aerothermal  prediction 


Figure  5.4.  Aerothermal  prediction  error  compared  to  midpoint  data 

To  focus  efforts  on  the  transient  aspects  of  the  aerothermal  problem,  the  connection  between 
the  stationary  and  transient  sections  of  the  Bayesian  network  is  supplemented  directly  by  the 
Glass  and  Hunt  pressure  datay^,  as  demonstrated  in  Figure  5.2.  This  removes  £PT  and  £yp from 

calibration  consideration  and  $  becomes  the  set  of  unknown  parameters  in  eERTi  and  £\  n  ,  which 
are  discrepancy  models  for  Eckert’s  reference  temperature  method  and  heat  transfer, 
respectively. 
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5.2  Global  and  Incremental  Discrepancy  Modeling  for  Time-Dependent  Problems 

In  a  related  study  by  Smarslok  et  al.  [23],  systematic  error  models  were  chosen  to  represent 
the  inadequacy  in  the  aerodynamic  pressure  and  heat  flux  predictions.  It  was  observed  that  the 
systematic  error  models  were  not  fully  sufficient  to  predict  the  data  with  reasonable  confidence 
across  the  dome.  Subsequently,  DeCarlo  et  al.  [24]  expanded  this  framework  to  include  quadratic 
error  models  with  respect  to  dome  slope  for  aerodynamic  pressure  and  heat  flux  predictions 
across  the  dome  to  capture  the  relationship  between  model  discrepancy  and  slope  of  the  dome. 
The  present  study  builds  upon  the  previous  framework  by  addressing  coupled  aerothermal  model 
discrepancies  through  time  across  the  dome. 

The  model  discrepancy  fonnulation  followed  in  this  study  is  based  on  the  Kennedy  and 
O'Hagan  (KOH)  calibration  framework  [21].  Relevant  sources  of  uncertainty  identified  in  the 
Bayesian  network,  including  data  uncertainty,  input  parameter  uncertainty,  and  input-specific 
model  discrepancy  have  representations  in  the  context  of  KOH  discrepancy,  shown  in  Eq.  (5.3). 

y(x,  t)  =  r](x,  t,  0)  +  s(x,  t,  </>)  +  ey  (5.3) 

where,  ;/  is  the  model  under  consideration,  </>  is  the  set  of  parameters  being  calibrated,  e  is  the 
model  discrepancy  function  evaluated  at  specific  input  vectors  x  at  time  t,  and  ev  ~  N( 0,ern)  is  the 
measurement  uncertainty  associated  with  that  set  of  data. 

The  transient,  coupled  aerothennal  solution  is  a  quasi-static  approximation  at  a  rate  of  At, 
which  lends  itself  to  two  ways  of  implementing  the  time-dependent  discrepancy  functions  in  the 
Bayesian  network.  The  first  strategy  involves  treating  model  discrepancy  globally  by  correcting 
the  time-dependent  quasi-static  prediction  after  the  set  of  n  evaluations  to  the  last  time  point  of 
interest  at  nAt.  Applying  the  discrepancy  function  in  this  way  offsets  the  prediction  by  s,  where 
the  fonn  can  be  estimated  a  priori  by  calculating  the  difference  between  the  n  uses  of  the  model 
and  data.  For  example,  the  discrepancy  between  a  linear  quasi-static  temperature  prediction  and 
Glass  and  Hunt  temperature  data  has  a  quadratic  trend  as  observed  in  Figure  5.4,  and  Eq.  (5.4) 
was  selected. 

£Ht  —  d0  +  d-it  -I-  d2t2  (5.4) 


The  discrepancy  in  Eckert’s  reference  temperature  follows  a  quadratic  trend  with  respect  to 
dome  slope,  as  demonstrated  by  DeCarlo  et  al.  [24],  In  addition  to  dependence  on  slope,  the 
discrepancy  in  heat  flux  must  also  account  for  a  functional  variation  through  time.  Thus,  a  linear 
relationship  with  time  is  assumed  for  £ert  in  Eq.  (5.5). 


dw 

£ert  —  co  +  cit  +  c2  +  c3 


(5.5) 


The  alternative  to  a  global  treatment  of  errors  is  by  incorporating  the  discrepancy  models  for 
Eckert’s  reference  temperature  and  linear  heat  transfer  in  Eqs.  (5.4)  and  (5.5)  incrementally  at 
each  time  step.  This  is  accomplished  by  correcting  the  prediction  each  time  one  of  the  coupled 
aerothennal  models  are  evaluated  and  recognizing  the  magnitude  of  errors  may  change 
throughout  the  time  domain.  Equation  (5.3)  corresponds  to  the  global  model  discrepancy  form.  A 
mathematical  representation  of  incremental  discrepancy  modeling  within  the  KOH  calibration 
framework  is  shown  in  Eq.  (5.6),  where  data  v,  at  each  time  step  At  is  explained  by  model 

41 

Approved  for  public  release;  distribution  unlimited. 


prediction  rjt  as  a  function  of  the  previous  model  output  77,-1,  model  discrepancy  e  and 
measurement  error  sr 


yi{x,t{)  =  +  £(*ftif0)  +  £y  (5.6) 

Figures  5.5  and  5.6  graphically  show  the  prediction  and  correction  for  global  and  incremental 
model  discrepancy  in  Eqs.  (5.3)  and  (5.6),  where  both  treatments  implement  the  time  and  slope- 
dependent  discrepancy  models  in  Eqs.  (5.4)  and  (5.5).  The  global  and  incremental  discrepancy 
strategies  are  used  to  correct  the  time-dependent  temperature  output.  For  a  global  model 
discrepancy  implementation,  the  model  output  is  a  full  time  history  up  until  nA t,  whereas 
incremental  discrepancy  treats  the  model  output  as  the  model  evaluation  at  every  At. 


Figure  5.5.  Correcting  Predictions  using  Global  Model  Discrepancy 


Figure  5.6.  Correcting  Predictions  using  Incremental  Model  Discrepancy 

The  incremental  approach  is  particularly  important  for  coupled,  multi-disciplinary  models, 
where  errors  are  being  propagated  between  models  at  each  time  step.  Figures  5.5  and  5.6  can  be 
re-imagined  for  the  coupled  aerothermal  models  where  Eckert’s  reference  temperature  is 
predicting  the  temperature  gradient,  while  linear  heat  transfer  uses  information  on  the  gradient 
and  predicts  the  next  temperature  in  time.  Correcting  coupled  models  incrementally  will  in  turn 
couple  their  model  discrepancies  and  exploit  their  relationship  through  the  Bayesian  network  to 
update  information  on  both. 

Further,  the  incremental  approach  results  in  separable  models  with  quantified  errors  and  can 
be  applied  individually  to  other  physical  systems  that  need  not  be  coupled.  This  is  in  contrast  to 
global  discrepancy  application,  where  the  calibration  corrects  the  cumulative  effect  of  errors  over 
time  as  a  consequence  of  using  both  models  indiscriminately.  Thus  the  coupled  aerothermal 
models  calibrated  using  global  discrepancy  must  be  taken  forward  as  a  unit. 
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Note,  the  incremental  and  global  discrepancy  functions  have  the  same  fonn,  however  the 
calibrated  coefficients  for  heat  transfer  discrepancy  (do,  d\,  di)  will  be  different  depending  on  the 
strategy  implemented.  Since  the  incremental  discrepancy  strategy  corrects  the  prediction  at  each 
time  step,  then  the  magnitude  of  the  discrepancy  at  a  given  time  t  is  less  than  that  of  the  global 
discrepancy.  Section  5.3  will  list  appropriate  prior  distributions  for  Bayesian  calibration  of  the 
model  discrepancy  parameters,  compare  the  incremental  and  global  error  application  approaches 
through  both  sensitivity  analyses  and  calibration,  and  select  an  error  treatment  to  apply  to  the  full 
dome  for  validation.  Section  5.3  will  list  appropriate  prior  distributions  for  Bayesian  calibration 
of  the  model  discrepancy  parameters. 

5.3  Prior  Uncertainty  and  Discrepancy  Implementation  Selection 

Prior  distributions  for  parameters  of  cert  and  cht  using  both  incremental  and  global  treatments 
of  model  discrepancy  are  shown  in  Table  5.1.  All  prior  distributions  are  uniform  distributions 
with  lower  and  upper  bounds  specified.  Prior  distribution  ranges  were  chosen  by  passing  the 
maximum  likelihood  values  of  the  other  parameters  through  the  transient  aerothermal  models 
and  comparing  nominal  predictions  to  the  Glass  and  Hunt  data. 

Table  5.1.  Prior  Distributions  for  Discrepancy  Model  Parameters 


Parameter 

Incremental 

Global 

Lower 

Bound 

Upper 

Bound 

Lower 

Bound 

Upper 

Bound 

Co  (W/cm2) 

-0.5 

0.3 

-0.5 

0.3 

Ci  (W/cm2f) 

-0.1 

0.1 

-0.1 

0.1 

Ci  (W/cm2) 

-100 

200 

-100 

200 

C3  (W/cm2) 

-1000 

1000 

-1000 

1000 

rfo(K) 

-0.1 

0.1 

-3 

3 

dA  (K It) 

-0.1 

0.1 

-4 

3 

d2  (K It2) 

-0.1 

0.1 

-2 

1 

Both  discrepancy  treatments  are  applied  through  the  dynamic  Bayesian  network  and  used  to 
calibrate  discrepancy  models  at  the  midpoint  of  the  dome.  The  calibration  domain  consists  of 
time-dependent  temperature  data  up  until  4  seconds  into  Runs  30,  31,  and  32  of  Glass  and  Hunt 
turbulent  flow  experiments.  The  remaining  second  of  temperature  data  is  used  for  decision¬ 
making.  For  this  comparison  only  the  dome  midpoint  temperature  data  (where  dw/dx  is  zero)  is 
considered;  the  quadratic  slope  tenns  have  no  effect  in  this  single  point  calibration. 

Figures  5.7  and  5.8  show  the  calibrated  incremental  and  global  discrepancy  models  for  heat 
transfer  through  time.  Note  that  both  models  have  similar  quadratic  trends,  however  incremental 
discrepancy  mean  predictions  are  on  average  two  orders  of  magnitude  less  than  those  of  the 
global  discrepancy  mean  prediction.  Further,  the  uncertainty  in  the  global  discrepancy  model  is 
growing,  whereas  the  uncertainty  in  the  incremental  discrepancy  model  is  relatively  constant. 
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Mean  discrepancy - 95%  Co nf.  bounds 

Figure  5.7.  Calibrated  global  discrepancy  model  through  time 


Mean  discrepancy - 95%  Co  nf.  bounds 

Figure  5.8.  Calibrated  global  discrepancy  model  through  time 


The  predictions  from  both  calibrated  discrepancy  models  at  the  dome  midpoint  decides  which 
should  be  used  for  the  remainder  of  the  dome  points.  The  calibrated  discrepancy  models  are 
propagated  forward  through  the  network  for  Runs  30,  3 1,  and  32  predictions. 


(a) 


Time(s) 

(b)  (c) 


Mean  Prediction - 95%  Conf.  Bounds  ♦  Dome  Midpoint  Data 


Figure  5.9.  Calibrated  aerothermal  model  predictions  using  global  discrepancy  application 
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Figure  5.10.  Calibrated  aerothermal  model  predictions  using  incremental  discrepancy 

application 

Figures  5.9  and  5.10  compare  the  predictions  from  calibrated  aerothermal  models  for  both 
discrepancy  treatments.  Note  that  the  calibration  domain  from  0  to  4  seconds  is  predominantly 
linear,  whereas  the  validation  temperature  data  deviates  significantly  from  linearity.  The  global 
model  discrepancy  demonstrates  a  linear  prediction  continuing  into  the  validation  domain. 
However,  the  incremental  discrepancy  model  is  able  to  capture  the  nonlinear  trend  in  the 
calibration  domain  and  extrapolate  it  to  the  validation  domain.  Confidence  bounds  in  Figures  5.9 
and  5.10  indicate  more  confidence  in  global  discrepancy  modeling  prediction;  however, 
particularly  in  Runs  30  and  32  there  is  deviation  between  the  prediction  and  data  in  the 
extrapolation.  Bayes  factor  is  used  to  more  formally  compare  the  validation  performance  of  the 
global  and  incremental  strategies.  Equation  (5.7)  shows  the  Bayes  factor  for  model  selection. 


Pr(yrl<fti,t;)  _  I  PrCyrlftx,  OftO d</> 

Pr(yrl02-ti)  /  Pr(yr|02,tj)7r2(02)d</> 


To  compute  the  Bayes  factor,  the  likelihood  ratio  of  two  competing  models  is  computed, 
which  incorporates  both  accuracy  and  precision  into  a  decision  metric.  In  Eq.  (5.7),  let  the 
likelihood  corresponding  to  incremental  discrepancy  model  be  defined  by  parameters  (/>\  and 
joint  probability  distribution  n\,  and  let  the  global  discrepancy  model  be  defined  by  parameters 
(f>2  and  joint  probability  distribution  tt2.  Therefore,  when  B(t,)  is  greater  than  1,  the  incremental 
discrepancy  prediction  supports  the  data  better  the  global  discrepancy.  The  integral  form  in  Eq. 
(5.7)  includes  the  likelihood  function  of  the  data  supporting  the  prediction  Pr(vjr/>)  and  the 
probability  density  function  (PDF)  of  the  both  hypotheses,  n\((j)\)  and  Bayes  factors 

through  the  validation  domain  are  shown  in  Figure  5.11. 
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Time  (s) 


Figure  5.11.  Likelihood  ratio  between  incremental  and  global  discrepancy  approaches 

through  time 


The  Bayes  factors  over  the  validation  domain  in  Figure  5.11  are  consistently  much  greater 
than  one,  favoring  the  incremental  approach  and  following  approximately  a  linear  trend  through 
time  on  a  semi-logarithmic  scale.  Averaged  time-dependent  prediction  standard  deviations  a 
(inverse  confidence)  in  both  approaches  every  0.2  seconds  are  listed  in  Table  5.2.  The 
uncertainty  from  both  discrepancy  models  grow  through  time.  The  Bayes  factors  are  so  large  due 
to  severe  bias  in  the  global  discrepancy  model  in  the  validation  domain. 


Table  5.2.  Averaged  time-dependent  temperature  uncertainty  and  likelihood  ratio 


Time  (s) 

^incremental 

(K) 

^global  (K) 

4.0 

0.47 

0.36 

4.2 

0.67 

0.43 

4.4 

0.91 

0.51 

4.6 

1.20 

0.60 

4.8 

1.55 

0.70 

5.0 

1.95 

0.81 

Figure  5.11  and  Table  5.2  indicate  the  incremental  approach  is  more  consistent  with  the  time- 
dependent  temperature  data  especially  farther  away  from  the  calibration  domain.  Recall,  the 
same  heat  flux  and  heat  transfer  discrepancy  models  were  used  in  both  cases  but  by  applying 
discrepancy  incrementally  we  have  greater  extrapolation  ability.  In  the  next  section,  the 
incremental  discrepancy  modeling  approach  is  applied  to  the  remaining  dome  centerline  points 
for  calibration. 

5.4  Calibration  with  Incremental  Discrepancy  for  Dome  Predictions 

Bayesian  calibration  of  the  model  discrepancy  coefficients  was  perfonned  using  the  dynamic 
Bayesian  network  with  the  time-dependent  aerothermal  data  from  Figure  5.2.  The  data  from  all 
three  runs  was  considered  to  calibrate  the  model  discrepancy  from  Eqs.  (5.4)  and  (5.5)  at  each  of 
the  eleven  points  along  centerline  of  the  Glass  and  Hunt  dome.  The  incremental  discrepancy 
treatment  was  chosen  for  this  part  of  the  investigation  based  on  the  Bayes  factor  comparison  in 
Figure  5.11  and  Table  5.2.  For  the  Bayesian  calibration,  the  posterior  distributions  for  the 
uncertainty  discrepancy  parameters  were  estimated  using  104  samples  from  the  Markov  Chain 
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Monte  Carlo  (MCMC)  algorithm  called  slice  sampling.  Figures  5.12-5.14  show  the  corrected 
predictions  using  the  incremental  discrepancy  approach  for  Runs  30,  31  and  32,  respectively.  For 
all  three  cases,  the  predictions  shown  at  1  and  3  seconds  are  in  the  calibration  domain,  whereas  5 
seconds  is  in  the  validation  domain. 


Figure  5.12.  Calibrated  Run  30  prediction  across  dome  at  1,  3,  and  5  seconds 


Mean  Prediction - 95%  Conf  Bounds  *  Dome  Data 


Figure  5.13.  Calibrated  Run  31  prediction  across  dome  at  1,  3,  and  5  seconds 


- Mean  Prediction - 95%  Cont  Bounds  *  Dome  Data 

Figure  5.14.  Calibrated  Run  32  prediction  across  dome  at  1,  3,  and  5  seconds 

Figures  5.12-5.14  show  small  model  discrepancy  in  the  calibration  time  domain  (t  =  Is  and 
3s)  and  growing  extrapolation  uncertainty  in  the  validation  domain  ( t  =  5s).  This  is  analogous  to 
Figure  5.10,  where  uncertainty  increased  when  predicting  beyond  4  seconds.  Table  5.3  lists  the 
average  uncertainty  in  the  calibrated  predictions  for  each  experiment.  Run  31  shows  the  most 
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prediction  uncertainty  at  5  seconds,  where  the  95  percent  confidence  bounds  are  larger  than  those 
in  Runs  30  and  32. 


Table  5.3.  Calibrated  temperature  prediction  uncertainty  at  t=l,  3,  and  5  seconds  for  Runs 

30,  31,  and  32 


Time  (s) 

<730  (K) 

<731  (K) 

<732  (K) 

1.0 

0.10 

0.11 

0.10 

3.0 

0.21 

0.25 

0.21 

5.0 

0.24 

0.33 

0.24 

The  next  section  will  use  the  model  reliability  metric  to  assess  the  confidence  in  the  corrected 
predictions  using  incremental  model  discrepancy. 

5.5  Confidence  Assessment  of  Calibrated  Aerothermal  Models  Using  the  Model 
Reliability  Metric 

To  assess  the  confidence  in  the  calibrated  discrepancy  models,  a  model  reliability  metric  was 
explored  as  a  measure  of  model  predictive  capability.  Equation  (5.8)  shows  the  reliability  r 
quantified  as  the  probability  that  the  difference  between  the  model  prediction  ym  and  data  yj 
being  within  a  tolerance  <5  [51].  Thus,  the  remaining  model  uncertainty  after  calibration  is  used  to 
assess  this  model  reliability. 


r  =  Pr[-S  <  yd  ~  Vm  <  8]  (5-8) 

The  probability  in  Eq.  (14)  may  be  estimated  through  Monte  Carlo  simulation,  or  efficient 
reliability  analysis  methods  such  as  the  first-order  or  second-order  reliability  methods 
(FORM/SORM).  The  model  reliability  is  straightforward  in  interpretation.  A  tolerance  of  S=  IK 
is  used  in  the  following  reliability  analysis. 
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Figure  5.15.  Reliability  of  calibrated  aerothermal  models  across  the  dome  in  the  validation 

domain  ( t  =  4-5s) 


Figure  5.15  shows  greater  confidence  in  the  calibrated  aerothennal  models  towards  the  middle 
of  the  dome.  This  suggests  that  the  points  near  the  edge  of  the  dome  exhibit  behavior  not 
captured  in  the  discrepancy  models,  which  was  observed  in  Figures  5.12  and  5.14.  Furthermore, 
the  model  reliability  metric  decreases  as  the  uncertainty  grows  through  time. 
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5.6  Summary  and  Conclusions 

This  research  is  part  of  a  series  of  investigations  on  quantifying  the  confidence  in  coupled 
aerothennal  elastic  models  for  hypersonic  aircraft  structures.  Bayesian  model  calibration 
methods  were  extended  for  a  coupled,  time-dependent  problem.  A  dynamic  Bayesian  network 
for  coupled  aerothennal  models  (including,  aerodynamic  pressure  and  heat  flux,  and  transient 
heat  transfer)  and  time-dependent  data  was  developed.  Bayesian  model  calibration  was 
performed  for  model  discrepancy  using  time-dependent  temperature  data  from  historic 
aerothennal  experiments  conducted  at  NASA's  8-foot  High-Temperature  Tunnel.  Two  ways  to 
implement  the  dynamic  model  errors  were  investigated,  and  it  was  determined  that  applying 
incremental  errors  through  time  supplied  the  model  with  greater  extrapolation  ability.  Finally,  the 
incremental  discrepancy  strategy  was  extended  to  calibrate  the  dynamic  aerothennal  error 
models  at  eleven  points  along  the  dome  centerline.  The  model  reliability  metric  for  model 
validation  was  used  to  assess  spatial  and  temporal  confidence  in  the  aerothennal  predictions 
across  the  dome.  Confidence  in  calibrated  aerothennal  discrepancy  models  increased  towards  the 
dome  midpoint  and  over  time  through  the  validation  domain. 
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6  INVESTIGATING  MODEL  UNCERTAINTY  IN  AEROELASTIC  RESPONSE  OF 
THIN  PANELS 

Recent  research  efforts  have  focused  on  the  development  of  computational  frameworks 
enabling  the  prediction  of  coupled  fhud-thermal-structural  response  [4-16].  McNamara  and 
Friedmann  (2011)  provided  a  comprehensive  review  of  different  solution  strategies  for 
calculating  the  response  of  a  hot  structure  in  a  hypersonic  flow  [4],  A  stochastic  collocation 
approach  was  used  by  Lamorte  et  al.  (2014)  for  propagating  uncertainty  in  aerothennoelastic 
analysis  [52].  A  subsequent  study  expanded  on  uncertainty  propagation  in  aerothennoelastic 
analysis  for  hypersonic  vehicles  with  emphasis  on  assessing  the  impact  of  aerothennoelastic 
deformation  on  aerodynamic  heating  [14].  Culler  and  McNamara  (2011)  identified  two-way 
coupling  between  structural  defonnation  and  aerodynamic  heating  as  an  important  consideration 
in  the  aerothennoelastic  modeling  of  a  panel  [11].  Skujins  (2013)  developed  a  reduced-order 
modeling  methodology  for  unsteady  aerodynamics  based  on  linear  convolution  with  a  nonlinear 
conection  factor  applicable  from  subsonic  to  hypersonic  flow  speeds  [53].  The  nonlinear 
conection  factor  is  computed  using  data  from  CFD  simulations  and  kriging.  These  efforts 
underscore  the  importance  of  understanding  the  uncertainty  in  a  coupled  aerothennoelastic 
model,  in  particular  model  uncertainty. 

To  reduce  the  computational  cost  of  long-duration  dynamic  response  simulations,  the  solution 
of  the  coupled  problem  can  be  obtained  with  reduced  order  models  (ROMs).  However,  reduced 
order  models  introduce  model  uncertainty  in  the  prediction  due  to  solution  approximation  and 
model-fonn  errors.  For  example,  piston  theory,  developed  by  Lighthill  (1953),  provides  a  simple 
point-wise  relation  between  the  aerodynamic  pressure  and  the  surface  motion  [54].  The  strengths 
of  this  model  are  its  simplicity  and  computational  efficiency;  on  the  other  hand,  the  accuracy  of 
piston  theory  greatly  diminishes  in  the  presence  of  significant  three-dimensional  flow  effects, 
combinations  of  high  Mach  number  and  surface  inclination,  and  viscous  effects  [55-57]. 

Uncertainty  exists  due  to  imperfect  knowledge  of  the  system  behavior,  physical  variability, 
model  order  reduction,  assumptions  and  approximations,  and  the  limited  experimental  data 
available  for  model  calibration  and  validation.  Liang  and  Mahadevan  (2011)  developed  a 
systematic  error  quantification  methodology,  which  distinguished  various  sources  of  uncertainty 
that  are  inherent  in  the  predictions  of  computational  models  [45].  These  sources  of  uncertainty 
were  grouped  in  three  categories:  (1)  numerical  errors  arising  from  the  model  inputs, 
discretization  errors  from  finite  element  analysis  (FEA)  mesh  discretization,  and  surrogate  model 
prediction  errors;  (2)  uncertainty  quantification  errors  due  to  the  finite  number  of  samples  taken 
in  the  uncertainty  propagation  analysis;  and  (3)  model-form  errors  that  arise  from  imperfect 
modeling  of  physics. 

The  quantification  of  model-form  errors  has  been  most  often  used  to  gain  understanding  on  the 
predictive  capability  of  computational  models.  Kennedy  and  O’Hagan  (2001)  proposed  a 
Bayesian  model  calibration  framework  to  account  for  various  sources  of  uncertainty  in  model 
predictions  and  the  calibration  of  uncertain  inputs  [21].  In  their  treatment  of  model  uncertainty, 
the  discrepancy  is  added  to  the  model  prediction  in  a  non-intrusive  fonn.  Moser  and  Oliver 
(2013)  developed  a  different  strategy,  where  the  model  inadequacy  is  introduced  at  the  source  of 
the  error  and  calibrated  with  the  uncertain  input  parameters  [58].  The  application  they  considered 
for  this  approach  was  the  calibration  of  turbulence  model  parameters  of  RANS  CFD  models. 

The  approach  developed  in  the  current  study  seeks  to  enhance  the  aerodynamic  pressure 
predictions  of  piston  theory  with  a  model-form  error  model.  This  is  achieved  through  quantifying 
the  model  inadequacy  in  the  local  slope  of  the  panel  induced  by  boundary  layer  displacement 
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effects.  This  approach  will  be  referred  to  as  effective  slope  model.  The  advantage  of  this 
approach  is  that  it  uses  the  physics  built  into  piston  theory  as  a  foundation,  and  incorporates 
information  from  higher-fidelity  models  (CFD)  to  quantify  the  model  uncertainty  to  improve 
predictions,  while  retaining  a  reduced  computational  cost.  This  follows  a  similar  philosophy  as  in 
the  multifidelity  optimization  community,  but  with  the  objective  of  response  prediction.  For 
example,  March  and  Willcox  (2012)  presented  a  provably  convergent  multifidelity  optimization 
algorithm  that  uses  radial  basis  function  interpolation  to  capture  the  error  between  high-fidelity 
and  low-fidelity  functions,  and  then  the  error  was  added  to  the  low-fidelity  function  to  create  a 
surrogate  model  of  the  high-fidelity  function  in  the  neighborhood  of  a  trust  region  [59]. 

The  impact  of  model  uncertainty  in  the  aerodynamic  pressure  acting  on  a  two-dimensional  flat 
panel  subjected  to  hypersonic  flow  will  be  explored  from  both  single-  and  cross-disciplinary 
perspectives.  From  the  cross-disciplinary  point  of  view,  the  nonlinear  flutter  response  of  the 
panel  will  be  computed  and  the  effect  of  model  uncertainty  in  the  aerodynamic  pressure  on  the 
limit-cycle  oscillation  amplitude  of  the  panel  will  be  analyzed. 

An  outline  of  this  section  is  as  follows.  Section  6.1  describes  the  two  components  of  the 
aeroelastic  computational  model:  1)  the  aerodynamic  pressure  model  (i.e.,  piston  theory),  and  2) 
the  structural  ROM  and  its  fonnulation.  Next,  discussion  on  the  development  and  verification  of 
the  piston  theory  with  effective  slope  is  presented  in  Section  6.2.  Construction  and  results  from 
the  verification  of  the  ROM  developed  for  the  structural  component  of  the  aeroelastic  system  are 
shown  in  Section  6.3.  The  forward  propagation  of  uncertainty  in  the  aeroelastic  model  response 
is  performed  in  Section  6.4  and  the  LCO  amplitude  and  frequency  results  from  piston  theory  and 
the  effective  slope  model  compared.  Finally,  a  sensitivity  analysis  is  performed  to  assess  how 
many  modes  should  be  used  in  the  selection  of  CFD  cases  for  the  construction  of  the  effective 
slope  approach. 

6.1  Aeroelastic  Model  Definitions 

Hypersonic  aircraft  structures  are  subjected  to  intense,  coupled,  fluid-thermal-structural 
loading  during  high-speed  flight.  The  aeroelastic  model  components  of  the  fluid-structure 
interaction  are  shown  in  Figure  6.1.  The  hypersonic  flow  acting  on  the  structure  results  in 
aerodynamic  pressure  on  the  wetted  surface  of  the  panel.  This  leads  to  elastic  deformation  of  the 
panel  into  the  flow  field,  resulting  in  feedback  on  the  flow. 


Figure  6.1.  Aeroelastic  coupling 

The  aeroelastic  model  is  used  to  investigate  the  impact  of  model  uncertainty  on  nonlinear 
panel  flutter.  Flutter  is  an  aeroelastic  instability  where  the  amplitude  of  vibration  of  a  structural 
component  in  a  flow  field  increases  without  a  bound.  In  the  case  of  a  panel,  nonlinear  membrane 
stretching  provides  a  stabilizing  effect  that  restrains  the  panel  motion  to  a  bounded  amplitude  for 
limit  cycle  oscillations  (LCO).  Panel  flutter  not  only  provides  an  extreme  response  scenario  for 
this  coupled  system,  but  is  a  design  constraint  of  aerospace  structures. 

Figure  6.1  is  a  schematic  of  the  solution  of  the  aeroelastic  problem.  The  panel  displacement 
and  velocity,  w(x,?)  andw(x,t) ,  are  computed  using  a  structural  ROM.  The  panel  deformation 
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serves  as  a  boundary  condition  to  the  flow  problem,  for  which  oblique  shock  relations  are  used  to 
compute  the  pressure  after  the  shock,  and  piston  theory  to  obtain  the  pressure  at  the  defonned 
surface  of  the  panel.  Note  that  piston  theory  assumes  that  the  flow  is  inviscid.  This  pressure  is 
then  used  as  a  right-hand-side  term  (i.e.,  loading)  in  the  solution  of  the  ROM  equations  of 
motion.  A  brief  description  of  piston  theory  and  the  structural  ROM  methodology  used  in  this 
work  are  presented  in  the  following  subsections. 


Structural  Response  Prediction 

w:  Nonlinear  ROM 


8w(x,t)  JL  , 

[dt 

dw{x,t)  ^dW,(x) 

8x  tt  dx  ‘  ^  ' 

M  ij  .  i  /)  //,  l  A-;0  //,  •  K™  //  // 

+  Vj  Vi  n„  =Ft 

Aerodynamic  Pressure  Prediction 

p3:  Oblique  Shock  Relations 
p4\  Piston  Theory 


Figure  6.2.  Aeroelastic  solution,  panel  slope  and  velocity  is  transferred  to  piston  theory  and 
aerodynamic  pressure  is  transferred  to  the  structural  solution 


6.1.1  Aerodynamic  Pressure 

Consider  a  panel  section  on  the  forebody  of  a  representative  hypersonic  vehicle  configuration, 
as  shown  in  Figure  2.1  [1 1].  As  the  vehicle  is  subjected  to  a  hypersonic  flow,  an  attached  oblique 
shock  is  created  at  the  forebody  leading  edge  (location  ‘1’).  The  surface  before  and  after  the 
panel  is  assumed  to  be  rigid;  therefore,  the  inviscid  flow  properties  at  locations  2  and  3  are  the 
same.  These  flow  properties  are  obtained  using  oblique  shock  relations. 

High-fidelity  numerical  simulations  of  the  nonlinear  structural  response  and  the  complex 
hypersonic  flow  environment  are  computationally  very  expensive.  An  uncertainty  quantification 
(UQ)  study  requires  several  samples  from  a  computational  model,  motivating  the  need  for 
computationally  efficient  tools.  Therefore,  piston  theory,  which  provides  a  simple  point-function 
relation  between  the  aerodynamic  pressure  and  the  motion  of  the  panel,  will  be  used  for 
computing  the  steady  and  unsteady  aerodynamic  loads  on  the  structure.  Piston  theory  has  been 
observed  to  provide  reasonably  accurate  pressure  predictions  for  sufficiently  large  Mach 
numbers  and  as  long  as  the  magnitude  of  the  normal  component  of  fluid  velocity  never  exceeds 
the  speed  of  sound  in  the  undisturbed  fluid  [9].  However,  the  accuracy  of  piston  theory  greatly 
diminishes  in  the  presence  of  significant  three-dimensional  flow  effects,  combinations  of  high 
Mach  number  and  surface  inclination,  and  viscous  effects  [18-21].  The  third-order  piston  theory 
expression  is  shown  in  Eq.  (2.5).  In  order  to  obtain  the  aerodynamic  pressure,  the  defonnation  of 
the  panel  needs  to  be  computed.  This  will  be  achieved  with  the  structural  ROMs  described  in  the 
next  subsection. 

6.1.2  Structural  Reduced  Order  Model  (ROM)  Formulation 

The  structural  ROMs  considered  in  this  study  are  based  on  a  representation  of  the  nonlinear 
geometric  response  in  terms  of  a  set  of  basis  functions, 
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(6.1) 


M 

u(0=IX(0  vin) 

n= 1 

where  u(?)  represents  the  vector  of  displacements  of  the  finite  element  degrees  of  freedom,  v|/('!) 
are  specified,  constant  basis  functions,  and  qn  (t  )  are  the  time  dependent  generalized  coordinates. 

The  ROM  procedure  described  here  is  achieved  in  the  undefonned  configuration  Q0  for 
which  the  field  equations  are  shown  in  Eq.  (6.2).  Note  that  summation  is  implied  over  repeated 
indices. 


~^r(Fy  Sjk )  +  Po  bi  =  Po (6.2) 

axk 

Where  S  is  the  second  Piola-Kirchhoff  stress  tensor,  p0  is  the  density  with  respect  to  the 

reference  configuration,  and  b°  is  the  vector  of  body  forces,  all  of  which  are  assumed  to  depend 
on  the  position  XeQ0,  [60,61].  Furthermore,  in  Eq.  (6.3),  F  denotes  the  deformation  gradient 
tensor  of  components 


F  =^L 

lJ  dX: 


=<s5+Al 

J  dX; 


(6.3) 


where  Sjj  is  the  Kronecker  delta  and  u  =  x  -  X  is  the  displacement  vector,  x  being  the  position 
vector  in  the  deformed  configuration.  In  the  present  formulation  the  material  is  assumed  to  be 
linear  elastic,  so  S  and  the  Green  strain  tensor  E  satisfy 

S,-C,jkIEkl  (6.4) 

where  C  is  a  fourth  order  elasticity  tensor,  which  is  a  function  of  the  undefonned  coordinates  X  . 

Carrying  on  with  the  formulation,  assume  next  the  displacement  field  uj  in  the  continuous 
structure  in  the  form 


M 

u,  (X,t)  =  Yjrin  (t)  U\n)  (X)  for  i  =  1,2,3  (6.5) 

n—l 

where  U\ " 1  (X)  are  specified,  constant  basis  functions  satisfying  the  boundary  conditions  also  in 

the  undefonned  configuration  and  this  is  the  continuous  space  equivalent  of  the  discrete,  finite 
element  model,  representation  of  Eq.  (6.1). 

By  introducing  Eq.  (6.5)  in  Eqs.  (6.2)-(6.4)  and  enforcing  the  condition  that  the  error  be 
orthogonal  to  the  basis  (i.e.,  Galerkin  approach),  a  set  of  nonlinear  ordinary  differential  equations 
for  the  generalized  coordinates  rjn  (t  )  can  be  obtained.  This  leads  [62]  to  the  reduced  order 
model  equations 
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+  D‘i'h  +  K"]  rh  +KS ^jVi+KfiVjV^p  =  Ft 


(6.6) 


where  a  linear  damping  term  Dtj  t)  j  has  been  added  to  collectively  represent  various  dissipation 
mechanisms.  Furthennore,  M..  denotes  the  elements  of  the  mass  matrix,  K]' 1 ,  are 

the  linear,  quadratic,  and  cubic  stiffness  coefficients  and  F.  are  the  modal  forces.  The  indirect 

evaluation  of  the  coefficients  in  Eq.  (6.6)  by  Hollkamp  and  Gordon  (2005)  from  a  finite  element 
model  is  used  here  [63]. 

With  the  aerodynamic  pressure  and  structural  response  models  defined,  the  next  section  will 
investigate  and  quantify  the  model-form  uncertainty  in  the  aerodynamic  pressure  predictions 
obtained  with  piston  theory. 

6.2  Model  Uncertainty  in  Aerodynamic  Pressure 

The  objective  of  mathematical  models  is  to  make  predictions  about  the  behavior  of  a  system. 
However,  no  model  is  perfect  due  to  various  sources  of  uncertainty,  so  the  predictions  will  not 
exactly  equal  the  true  value  of  the  process  it  is  intended  to  represent.  Model  uncertainty  (also 
referred  to  as  model  inadequacy)  can  be  defined  as  the  difference  between  the  true  mean  of  the 
process  and  the  output  of  the  model  at  specified  inputs  [21].  Model  uncertainty  S  can  be 
introduced  in  two  general  forms.  First,  by  adding  a  non-intrusive,  external  model  discrepancy 
tenn,  as  shown  in  Eq.  (6.7). 


y  =  m{x,(t>)  +  5  +  e  (6.7) 

Where  jc  are  the  model  inputs,  (f)  are  the  uncertain  model  parameters,  and  s  represents  the 
measurement  uncertainty.  A  separate  model  discrepancy  term  is  a  commonly-used  and  flexible 
approach  to  quantifying  the  model  fonn  uncertainty  in  m  [21,45].  However,  it  may  be  desirable 
to  capture  the  nature  of  the  model  inadequacy  at  the  source  of  the  error  within  the  model  itself,  as 
shown  in  Eq.  (6.8)  [58].  This  can  help  preserve  the  underlying  physics  and  original  model  form, 
while  still  providing  an  improved,  uncertainty-quantified  prediction. 

y  =  m(x,<j>,S)  +  s  (6.8) 

The  external  and  internal  model  discrepancy  approaches  are  discussed  in  the  following 
subsections  in  the  context  of  model  uncertainty  in  aerodynamic  pressure  predictions  using  piston 
theory. 

6.2.1  External  Discrepancy  Model 

Under  realistic  flow  conditions,  model-form  error  exists  in  piston  theory  due  to  unmodeled 
flow  viscosity  and  significant  three-dimensionality  of  the  flow  [55-57].  Using  data  available 
from  tests  conducted  by  Glass  and  Hunt  (1986)  on  spherical  dome  protuberances  subjected  to 
hypersonic  flow  [44],  Culler  and  McNamara  (2010)  [10]  showed  peak  and  average  pressure 
discrepancies  in  piston  theory  of  54.5%  and  18.6%,  respectively.  In  a  different  study,  DeCarlo  et 
al.  (2013)  identified  and  calibrated  through  a  Bayesian  network  a  model-form  error  function  with 
respect  to  the  slope  of  the  surface  for  selected  domes  from  the  Glass  and  Hunt  tests  [24].  The 
fonn  of  the  model  (i.e.,  quadratic  polynomial)  was  inferred  from  an  initial  interrogation  of  the 
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errors  between  3ld-order  piston  theory  predictions  and  experimental  observations.  Perez  et  al. 
(2013)  went  on  to  incorporate  the  calibrated  discrepancy  model  terms  into  piston  theory,  and 
used  the  resulting  model  for  the  prediction  of  flutter  and  LCO  amplitude  of  a  panel  in  the  post¬ 
flutter  regime  [32],  In  addition,  a  preliminary  verification  study  was  conducted  to  assess  the 
quality  of  the  calibrated  discrepancy  model  with  a  larger  range  of  slopes.  The  data  considered  for 
verification  were  CFD  results  published  by  Nydick  et  al.  (1995)  [64],  who  compared  piston 
theory  to  full  Navier-Stokes  CFD  simulations  using  a  two-dimensional  panel  with  a  prescribed 
wall  motion.  Third-order  piston  theory  with  the  calibrated  discrepancy  model  from  [24]  showed 
significant  improvement  compared  to  the  original  piston  theory  predictions. 

Verification  of  the  discrepancy  model  is  revisited  in  this  current  research  by  comparing 
pressure  predictions  to  higher-fidelity  data  obtained  from  a  kriging  surrogate  constructed  from  a 
set  of  Navier-Stokes  (i.e.  RANS)  CFD  simulations  of  flow  on  a  2-D  panel  [55].  The  structural 
deformations  in  the  construction  of  the  kriging  surrogate  were  parameterized  in  terms  of  the  first 
six  structural  mode  shapes.  The  flow  and  structural  parameters  used  in  the  construction  of  the 
surrogate  are  shown  in  Table  6.1.  The  panel  is  assumed  to  be  simply-supported  and  located  on  a 
rigid  wedge,  as  shown  in  Figure  2.1.  Structural  deformations  of  the  panel  proportional  to  the 
first,  second,  and  third  mode  shapes  were  enforced  for  different  modal  amplitudes  and  the 
pressure  coefficient  computed  with  piston  theory  (PT),  piston  theory  with  the  discrepancy  model 
(PTs),  and  the  surrogate  model  (Suit).  The  first  three  mode  shapes  are  shown  in  Figure  6.3.  The 
absolute  error  between  Cspurr  and  CPT  as  a  function  of  slope,  and  the  values  of  CSurr ,  CPT  ,  and 

CPpTs  as  a  function  of  the  location  along  the  streamwise  direction  are  shown  in  Figures  6.4-6. 6  for 
different  modal  amplitudes.  These  results  show  that  the  difference  between  C1’1  and  Cpun  can  be 

modeled  as  a  function  of  slope  for  a  defonnation  proportional  to  mode  1  (i.e.,  similar  to  the  Glass 
and  Hunt  domes).  However,  for  higher-order  modes  (e.g.,  modes  2  and  3)  two  or  more  Cp  values 
are  possible  for  a  given  slope;  therefore  slope  alone  cannot  fully  describe  the  pressure  trends  for 
these  structural  configurations.  In  all  cases,  the  pressure  predicted  by  the  surrogate  model  shows 
a  sharp  drop  near  the  leading  edge  of  the  panel,  which  is  not  captured  by  the  piston  theory 
models.  Due  to  the  inviscid  flow  assumption  used  in  piston  theory,  the  non-zero  slope  at  the  edge 
of  the  panel  becomes  a  sharp  discontinuity  for  the  flow.  Nevertheless,  the  boundary  layer  leads 
to  smooth  flow  from  the  flat  rigid  surface  of  the  deformed  panel. 

Table  6.1.  Fluid  and  structural  parameters  for  initial  verification  study 


Parameter 

Value 

Unit 

Altitude 

30 

km 

Freebody  Surface 

s 

deg 

Inclination 

Freestream  Mach 

8 

N/A 

Number 

Panel  Length 

1.5 

m 

Panel  Thickness 

5 

mm 

Twall 

300 

K 
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Figure  6.3.  Transverse  displacement  component  of  the  first,  second,  and  third  mode  shapes 
of  a  2-D  panel.  Each  mode  is  scaled  by  a  modal  amplitude  a,  for  i  =  1,  2,  and  3. 


(a)  (b) 

Figure  6.4.  Pressure  coefficient  for  panel  displacement  proportional  to  mode  1  with  peak 
displacement  equal  to  10  panel  thicknesses  at  M\  =  8:  (a)  pressure  coefficient  vs.  location 
along  the  streamwise  direction,  (b)  absolute  error  between  Cf"T  and  CPT  vs.  slope. 
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(a)  (b) 

Figure  6.5.  Pressure  coefficient  for  panel  displacement  proportional  to  mode  2  with  peak 
displacement  equal  to  3  panel  thicknesses  at  M\  =  8:  (a)  pressure  coefficient  vs.  location 
along  the  streamwise  direction,  (b)  absolute  error  between  Cspurr  and  CppT  vs.  slope. 


(a)  (b) 

Figure  6.6.  Pressure  coefficient  for  panel  displacement  proportional  to  mode  3  with  peak 
displacement  equal  to  1.5  panel  thicknesses  at  Mi  =  8:  (a)  pressure  coefficient  vs.  location 
along  the  streamwise  direction,  (b)  absolute  error  between  Cspurr  and  CpT  vs.  slope. 

The  results  presented  in  this  section  showed  that  the  discrepancy  in  the  pressure  is  not  only  a 
function  of  the  slope,  but  also  a  function  of  the  location  along  the  streamwise  direction  on  the 
panel.  This  observation  will  be  used  in  the  next  subsection  in  the  development  of  an  error  model 
that  incorporates  information  from  the  viscous  effects  not  modeled  in  piston  theory.  This  new 
approach  differs  from  the  discrepancy  model  presented  in  this  section  in  that  the  model-fonn 
error  is  introduced  at  the  source,  rather  than  externally. 

6.2.2  Internal  Error  Model  using  Effective  Slope 

Boundary  layer  displacement  effects  have  been  identified  as  a  potential  issue  with  piston 
theory  predictions  in  the  literature  [55-57].  Hypersonic  flows  are  characterized  by  relatively 
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thick  boundary  layers  [64]  which  displace  the  outer  inviscid  flow  and  cause  the  body  shape  to 
appear  much  thicker.  Accounting  for  the  boundary  layer  displacement  effects  can  be 
accomplished  by  the  introduction  of  an  effective  slope.  Then,  the  steady  coefficient  of  pressure 
predicted  from  piston  theory  is  given  by  Eq.  (6.9) 


c;rw= 


M. 


dw 


dx 


« +r±V 


dw 


eff 


v  dx  j 


—Ml 

12 


dw 


eff 


v  dx  j 


(6.9) 


where  dweff/dx  is  the  effective  slope.  The  effective  slope  is  detennined  using  the  approach 

proposed  by  McNamara  (2005)  [65],  where  the  coefficient  of  pressure  computed  from  a  CFD 
solution  of  the  steady  Navier-Stokes  equations  is  equated  with  Eq.  (6.9)  at  every  location  x 
desired 


c;rM-cf(*)= 0 


(6.10) 


where  C^s  is  the  pressure  coefficients  obtained  from  the  CFD  solution.  The  use  of  steady  CFD 

data  is  based  on  the  assumption  that  the  steady  flow  components  are  dominant  in  hypersonic 
flows  [56]. 

The  concept  of  effective  slope  leads  to  the  notion  that  flow  viscous  effects  introduce 
uncertainty  in  the  local  slope  of  the  panel.  The  quantification  of  this  uncertainty  is  propagated 
forward  to  determine  the  model-form  uncertainty  in  piston  theory.  This  is  performed  using 
infonnation  on  geometric  configurations  obtained  from  the  structural  dynamics  (i.e.,  the  mode 
shapes  of  the  structure).  The  uncertainty  in  the  local  slope  of  the  panel  is  represented  in  Eq. 
(6.1 1)  for  a  defonnation  proportional  to  the  structural  mode  shape  i 

dwiffiX’dnP^M^)  _  d(p\x 
dx  dx 


) 


(6.11) 


where  dweff  jdx  is  the  effective  slope  which  takes  into  account  the  boundary  layer  displacement, 
dcpjdx  is  the  physical  slope  of  mode  i,  tp  is  the  corresponding  weight  or  modal  amplitude,  and 
S  is  the  error  model.  The  error  model  is  a  function  of  space  x,  the  modal  amplitude  //„  and 
freestream  flow  properties  p-,.  and  M,..  A  general  structural  defonnation  can  be  parameterized  in 
terms  of  a  set  of  N  mode  shapes.  For  this  case,  the  effective  slope  is  represented  as 


dweff(x,ril,...,rii,px,Mx) 

dx 


n  ( 

s 


/■=!  V 


d<Pj{x) 

dx 


\ 

Vi  +  <5;  (•*,  r/i ,  Poo ’ ) 

) 


(6.12) 


The  error  model  S,  contains  zeroth-order  effects  that  come  from  boundary  layer  information  of 

the  undeformed  panel.  To  avoid  overprediction  of  the  effective  slope,  the  zeroth-order  effects  in 
Eq.  (6.12)  are  only  accounted  for  in  the  mode  1  enor  model. 

As  seen  in  Eqs.  (6.13)  and  (6.14),  the  error  models  is  a  function  of  several  variables.  In  order  to 
model  these  dependencies,  8t  is  obtained  by  fitting  the  error  between  the  effective  and  physical 
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slopes  using  kriging.  It  is  important  to  highlight  that  the  effective  slope  model  needs  a  source  of 
high-fidelity  data,  but  it  is  not  restricted  to  computational  data;  wind  tunnel  test  data  could  also 
be  used. 

A  kriging  model  is  made  of  two  components,  a  regression  function  fix)  is  constructed  using 
the  data  and  a  Gaussian  process  Z(x )  is  built  from  the  residuals  [66].  This  is  expressed 
mathematically  as 


Y(x)  =  f(x)  +  Z(x)  (6.13) 

where  the  Gaussian  process  Z(x)  has  a  mean  of  zero,  variance  a  ,  and  a  correlation  matrix  W. 

It  is  worth  noting  that  the  error  model  could  also  be  constructed  for  a  combination  of  modes 
instead  of  individually  for  every  mode  in  the  basis  as  follows 


dweff(x,Jil,...,Jii,pta,Mx)  d<p,(x) 


dx 


I 

1=1 


dx 


■7 , 


+  8(x,Tjl,...,Tji,paa,Mai) 


(6.14) 


The  fonn  of  the  error  model  shown  in  Eq.  (6.14)  will  be  used,  as  it  simplifies  the  construction  of 
the  kriging  models  by  reducing  the  size  of  the  input  space. 

It  is  assumed  that  the  primary  features  of  the  flow  are  captured  using  a  steady-state  analysis 
[55].  The  unsteady  effects  due  to  the  surface  motion  of  a  vibrating  panel  are  accounted  for  by 
incorporating  the  time-dependent  terms  of  3ld-order  piston  theory. 

The  error  model  Si  was  constructed  using  the  Fortran  Kriging  (ForK)  Library  with  a  2nd-order 
polynomial  and  the  Matem  correlation  function  with  v  =  3/2.  The  hyperparameters  of  the  kriging 
model  were  determined  using  a  pattern  search  algorithm.  Error  models  for  the  first  four  structural 
modes  were  constructed  for  the  parameter  space  shown  in  Table  6.2.  The  selection  of  the  range 
for  each  modal  amplitude  was  based  on  the  expected  structural  deformations  from  a  preliminary 
aeroelastic  analysis  of  the  panel  undergoing  LCO  with  the  aerodynamic  pressure  computed  with 
piston  theory. 

Table  6.2.  Discrepancy  model  parameter  space. 


Parameter 

Values 

Unit 

Altitude 

29,30,31 

km 

Freestream  Mach 
Number 

5,7,  9,  12 

N/A 

Twall 

300 

K 

a\ 

±4.0,  3.0,  2.0, 
1.0 

Non-dim 

<72 

±4.0,  3.0,  2.0, 
1.0 

Non-dim 

a?, 

±3.0,  2.0,  1.0, 
0.5 

Non-dim 

<74 

±1.5,  1.0,  0.5 

Non-dim 

The  high-fidelity  data  used  for  the  construction  and  verification  of  the  effective  slope  error 
model  were  obtained  from  solutions  of  the  Navier-Stokes  equations  using  the  NASA  Langley 
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CFL3D  code.  This  code  uses  an  implicit,  finite-volume  algorithm  based  on  upwind-biased  spatial 
differencing  to  solve  the  Reynolds  Averaged  Navier-Stokes  (RANS)  equations.  No  real  gas 
effects  were  included  in  the  analysis  since  CFL3D  is  an  ideal  gas  code.  The  Menter  k-co  SST 
turbulence  model  is  used  here  for  closure  of  the  RANS  equations  [39].  The  computational  mesh 
used  in  this  study  is  the  one  built  by  Crowell  et  al.  (2011)  with  a  total  of  34,600  cells  with  a 
maximum  v+  value  at  the  wall  of  0.65  [55].  Also,  the  upstream  surface  was  modeled  such  that 
the  transition  to  turbulent  flow  occurs  one  meter  upstream  of  the  leading  edge  of  the  panel.  A 
clamped-clamped  2-D  panel  is  considered  in  this  study  as  these  boundary  conditions  are  more 
representative  of  what  is  found  in  real  aircraft  structures. 

The  first  four  mode  shapes  of  the  2-D  clamped-clamped  panel  are  shown  in  Figure  6.7;  note 
that  the  zero  slope  at  the  leading  and  trailing  edges  of  the  panel  eliminate  the  discontinuity  seen 
by  the  flow  that  was  discussed  for  the  simply-supported  panel.  Shown  in  Figure  6.8  is  the 
pressure  coefficient  along  the  length  of  the  panel  for  rigid  panel  shape  obtained  from  the 
combination  of  the  first  four  structural  modes  with  the  following  modal  amplitudes:  a\  =  +3,  cii  = 
-3.25,  <23  =  +1.25,  and  <74  =  -0.15.  The  results  were  obtained  at  a  freestream  Mach  number  of  M\ 
=  10  and  an  altitude  of  30km,  and  computed  with  CFD  (NS),  piston  theory  (PT)  and  piston 
theory  with  effective  slope  ( PTeff ).  The  PTeff  results  show  a  close  agreement  with  the  NS  data. 

Based  on  the  Lx  error  metric,  the  error  between  CPTeff  and  CNS  is  of  4.7%  versus  30%  between 

Cp  and  C  .  The  remaining  model  uncertainty  in  the  C  ^prediction  is  likely  to  come  from 

the  limitations  of  kriging  and  the  use  of  a  finite  number  of  training  data. 

The  verification  under  unsteady  conditions  was  perfonned  by  enforcing  the  panel  motion  as 
follows 


w(x,t)  =  <2;  sin(W)<^.  (x)  (6.15) 

where  a,  is  the  amplitude  of  mode  (p,(x)  and  co  is  the  frequency  of  vibration  in  rad/sec.  Shown  in 
Figures  6.8-6.12  are  the  Generalized  Aerodynamic  Forces  (GAFs)  obtained  from  CFD 

( GAFt'NS) ),  piston  theory  ( GAF{IT] ),  and  piston  theory  with  effective  slope  ( GAF ^PTeff^),  From 
Eq.  (6.16),  L  and  c  are  the  length  and  width  of  the  panel  for  displacements  proportional  to  modes 
1,  2,  3,  and  4,  an  oscillation  frequency  of  100Hz,  and  the  same  freestream  flow  properties 
considered  in  the  verification  of  the  steady-state  analysis.  The  modal  amplitudes  of  the  panel 
motion  used  were  ai  =  +3.0,  <22  =  +2.0,  as  =  +1.0,  and  <24  =  +0.5.  Table  6.3  shows  results  that 

quantitatively  compare  GAF{NS')  to  GA F' P1 1  and  GAF^PTeff\  With  maximum  errors 

approximately  an  order  of  magnitude  less  than  GAF(PT\  the  GAF results  show  closer 
agreement  with  GAF1  VVi . 

Recall  that  the  intent  of  this  research  is  to  identify  missing  physics  within  computational 
models,  and  develop  an  internal  error  model  to  quantify  the  model-form  uncertainty.  A  kriging 
surrogate  was  selected  to  model  the  internal  error  through  effective  slope.  This  is  in  contrast  to 
external,  non-intrusive  approaches,  as  those  represented  by  Eq.  (6.8),  where  the  model-form 
uncertainty  is  quantified  directly  by  comparing  the  model  predictions  against  experimental 
observations.  In  the  case  of  the  PT  model,  viscous  effects  in  the  form  of  boundary  layer 
displacement  thickness  were  identified  as  important  information  that  is  not  included  in  the 
model.  The  objective  of  the  PTeff  model  was  to  incorporate  this  missing  infonnation  in  PT.  The 
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quantification  of  model-form  uncertainty  with  the  PTejf  model  corroborates  the  importance  of 
viscous  effects  in  hypersonic  flows. 


L 

GAFj  (t)  =  cjft  (x)p4  (x,t)  dx  for  i=l,  (6.16) 

0 


Table  6.3.  Difference  in  GAFs  for  the  Navier-Stokes  solution  compared  to  piston  theory 
and  piston  theory  with  the  effective  slope  model. 


Max  Diff. 

Max  Diff. 

Max  Diff. 

Max  Diff. 

Method 

GAFs 

GAFs 

GAFs 

GAFs 

Mode  1 

Mode  2 

Mode  3 

Mode  4 

(%) 

(%) 

(%) 

(%) 

gaf(pt) 

120 

100.3 

107.6 

100.2 

GAF 

10.6 

9.2 

8.6 

7.6 

x/L 

Figure  6.7.  Transverse  displacement  component  of  the  first,  second,  third,  and  fourth  mode 
shapes  of  a  2-D  clamped-clamped  panel.  Each  mode  is  scaled  by  a  modal  amplitude  «,•  for 

i=l,  2,  3,  and  4 
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0.04 


Figure  6.8.  Navier-Stokes,  piston  theory,  and  piston  theory  with  effective  slope  model 
pressure  coefficient  for  panel  displacement  proportional  to  a  combination  of  modes  1,  2,  3, 
and  4  at  Mi  =  10,  with  modal  amplitudes  a?  =  +3.0,  «2  =  -3.25,  a 3  =  +1.25,  CI4  =  -0.15 


Figure  6.9.  Navier-Stokes,  piston  theory,  and  piston  theory  with  effective  slope  Generalized 
Aerodynamic  Forces  for  enforced  panel  motion  w(x,t)  =  ai  sin(cot)(p(x)  proportional  to  mode 

1  and  frequency  equal  to  100Hz  at  Mi  =  10 
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Figure  6.10.  Navier-Stokes,  piston  theory,  and  piston  theory  with  effective  slope 
Generalized  Aerodynamic  Forces  for  enforced  panel  motion  w(x,t)  =  a2  sin(cot)(p(x) 
proportional  to  mode  2  and  frequency  equal  to  100Hz  at  Mi  =  10 


Non-Dimensional  Time  (Cycles) 

Figure  6.11.  Navier-Stokes,  piston  theory,  and  piston  theory  with  effective  slope 
Generalized  Aerodynamic  Forces  for  enforced  panel  motion  w(x,t)  =  (13  sin((ot)ip(x) 
proportional  to  mode  3  and  frequency  equal  to  100Hz  at  Mi  =  10 
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Figure  6.12.  Navier-Stokes,  piston  theory,  and  piston  theory  with  effective  slope 
Generalized  Aerodynamic  Forces  for  enforced  panel  motion  w(x,t)  =  (14  sin(cot)<p(x) 
proportional  to  mode  4  and  frequency  equal  to  100Hz  at  Mi  =  10 

The  uncertainty  quantified  in  this  model  will  be  propagated  to  the  aeroelastic  quantity  of 
interest  (i.e.,  limit  cycle  oscillations)  in  Section  6.4.  However,  the  next  step  is  to  generate  and 
verify  the  other  component  to  the  coupled  aeroelastic  model,  namely  the  structural  ROM. 

6.3  Structural  Reduced-Order  Model  (ROM)  Construction  and  Verification 

A  structural  ROM  introduces  uncertainty  in  the  analysis  in  the  fonn  of  solution  approximation 
errors,  which  is  related  to  the  number  of  basis  functions  used  to  construct  the  ROM.  As  described 
in  Section  6.1,  the  structural  ROMs  used  in  this  study  are  built  from  a  finite  element  analysis 
(FEA)  model;  therefore,  its  predictions  can  only  be  as  good  as  the  predictions  from  the  FEA 
model  that  was  used  to  construct  it. 

A  ROM  built  with  6  linear  modes  was  constructed  for  an  isotropic  2-D  panel  clamped  along 
the  sides  perpendicular  to  the  direction  of  the  flow.  The  geometric,  material  properties  of  the 
panel,  and  flow  conditions  are  shown  in  Table  6.4.  The  ROM  was  constructed  with  an  in-house 
FEA  beam  model  based  on  a  co-rotational  formulation  capable  of  analyzing  problems  with  large 
rotations  and  small  strains.  The  FEA  model  was  built  using  40  beam  elements,  a  total  of  123 
degrees-of-freedom.  A  mesh  convergence  analysis  was  performed  with  respect  to  the  first  six 
natural  frequencies  of  the  panel  and  the  nonlinear  stiffness  coefficient  Ann. 

Table  6.4.  Aeroelastic  model  parameters 


Parameter 

Value 

Unit 

Altitude 

30 

km 

Freebody  Surface  Inclination 

5 

deg 

Freestream  Mach  Number 

5-12 

N/A 

Panel  Length 

1.5 

m 

Panel  Thickness 

2 

mm 

Modulus  of  Elasticity 

113 

GPa 

Density 

4539 

kg/m3 

64 

Approved  for  public  release;  distribution  unlimited. 


The  ROM  and  the  FEA  model  were  coupled  with  PT  to  predict  the  aeroelastic  response  of  the 
panel.  The  quantities  of  interest  are:  the  LCO  amplitude  and  frequency.  These  quantities  are 
critical  for  aircraft  structural  design,  since  they  drive  the  fatigue  life  of  the  panel.  The  LCO 
amplitude  is  the  maximum  displacement  for  every  period  of  oscillation  and  it  occurs  at  the  14- 
point  along  the  length  of  the  panel. 

Shown  in  Figure  6.13  is  a  comparison  of  the  LCO  amplitude  as  a  function  of  the  Mach 
number.  The  maximum  error  of  the  6-mode  ROM  relative  to  the  FEA  results  occurs  at  a  Mach 
number  of  12  and  is  equal  to  approximately  2%.  Shown  in  Figure  15  is  the  LCO  frequency  as  a 
function  of  the  Mach  number.  The  maximum  relative  error  in  the  LCO  frequency  between  the  6- 
mode  ROM  and  the  FEA  is  equal  to  10%  which  corresponds  to  a  difference  of  approximately 
2Hz.  As  the  LCO  amplitude  grows,  oscillations  about  the  deformed  configuration  become  stiffer, 
this  is  seen  by  the  increasing  LCO  frequency  with  respect  to  the  Mach  number.  The 
quantification  of  the  solution  approximation  error  introduced  by  the  6-mode  ROM  suggests  that 
more  modes  would  be  needed  for  Mach  numbers  larger  than  12. 


Mach  Number 


Figure  6.13.  LCO  amplitude  at  panel  three-quarter  point  as  a  function  of  Mach  number 
for  air  properties  calculated  at  an  altitude  of  30km,  FEA,  6-Mode  ROM,  and  4-Mode  ROM 

coupled  with  3,d-order  piston  theory 
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Figure  6.14.  LCO  frequency  as  a  function  of  Mach  number  for  air  properties  calculated  at 
an  altitude  of  30km,  FEA,  6-Mode  ROM,  and  4-Mode  ROM  coupled  with  3' ''-order  piston 

theory 

The  next  section  investigates  the  effect  of  model-form  uncertainty  in  aerodynamic  pressure 
predictions  on  the  LCO  amplitude  and  frequency  of  a  panel. 

6.4  Uncertainty  Propagation  to  Coupled  Aeroelastic  Model  Response 

The  internal  error  from  Section  6.2  captured  the  effect  of  model-form  uncertainty  through  the 
effective  slope  used  in  aerodynamic  pressure  predictions  from  piston  theory.  The  PTeff  model  is 
used  here  to  investigate  the  effect  of  model-form  uncertainty  from  aereodynamic  pressure 
predictions  on  the  coupled  aeroelastic  response  of  the  2-D  panel  described  in  the  previous 
section.  The  structural  solution  was  obtained  using  the  6-mode  ROM;  however,  the  PTeff  model 
only  contains  information  from  the  first  four  dominant  modes,  as  the  modal  amplitudes  of  the 
last  two  modes  are  relatively  smaller.  The  validity  of  this  assumption  will  be  investigated 
through  the  sensitivity  analysis  of  the  four  modes  at  the  end  of  this  section. 

Variability  in  the  inputs  was  modeled  with  statistical  distributions  for  the  freestream  pressure  p\ 
and  the  modulus  of  elasticity  E.  Normal  distributions  were  used  with  means  equal  to  their 
nominal  values  shown  in  Table  6.4  and  1%  coefficients  of  variation.  Prior  knowledge  on  the 
model-form  uncertainty  in  PT  predictions  was  incorporated  based  on  previous  reports  [9]  which 
indicated  that  PT  predictions  are  expected  to  be  accurate  within  [-19%,  0%]. 

Figures  6.15  and  6.16  show  the  mean  and  95%  confidence  bounds  of  the  LCO  amplitude  and 
frequency  as  a  function  of  the  freestream  Mach  number  from  1,000  Monte  Carlo  samples.  The 
effects  of  model-form  uncertainty  in  aerodynamic  pressure,  quantified  with  the  PTeff  model, 
become  more  pronounced  for  increasing  Mach  numbers.  At  a  Mach  number  of  12,  the  mean 
values  of  the  LCO  amplitude  and  frequency  are  44%  and  13%  lower,  respectively,  when  model- 
form  uncertainty  from  viscous  effects  in  the  aerodynamic  pressure  prediction  is  incorporated. 
Increasing  discrepancies  in  the  LCO  amplitude  and  frequency  imply  larger  uncertainty  in  the 
stress  amplitude  and  number  of  cycles  that  the  structure  is  subjected  to.  Therefore,  future  work 
will  consider  the  analysis  of  the  impact  of  model  uncertianty  on  the  fatigue  life  of  the  panel. 
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Figure  6.15.  LCO  amplitude  at  panel  three-quarter  point  as  a  function  of  Mach  number 
for  nominal  air  properties  calculated  at  an  altitude  of  30km,  6-Mode  ROM  coupled  with 
3ld-order  piston  theory  and  3ld-order  piston  theory  with  effective  slope  model 


Figure  6.16.  LCO  frequency  as  a  function  of  Mach  number  for  nominal  air  properties 
calculated  at  an  altitude  of  30km,  6-Mode  ROM  coupled  with  3ld-order  piston  theory  and 
3l  d-order  piston  theory  with  effective  slope  model 

Next,  the  sensitivity  of  the  LCO  amplitude  and  frequency  to  the  aerodynamic  pressure 
contributions  from  the  first  four  modes  in  the  basis  is  assessed.  To  this  end,  distributions  of  the 
PTeff  model  predictions  were  constructed  by  adding  a  ±1%  error  to  the  PTejj  computations  for 
modes  1-4.  The  main  and  total  effect  indeces  were  computed  using  the  approach  described  by 
Saltelli  et  al.  (2010)  which  leads  to  a  reduction  of  the  computational  cost  from  N  to  N(k  +  2), 
where  N  represents  the  number  of  samples  and  k  is  the  number  of  uncertain  input  factors  [67]. 
The  exploration  of  the  input  space  is  performed  using  Sobol’  quasi-random  sequences  [68,69]  as 
suggested  by  Saltelli  [67].  The  results,  shown  in  Figures  6.17  and  6.18,  indicate  that  the  error 
model  of  mode  2  plays  a  very  important  role.  The  interaction  effects,  represented  by  the  total 
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effect  index,  are  stronger  than  the  first-order  effects  for  modes  1,3,  and  4.  These  results  suggest 
that  the  error  models  of  the  asymmetric  modes  (i.e.  modes  2  and  4)  play  a  more  important  role 
than  the  symmetric  modes  (i.e.  modes  1  and  3).  This  is  not  surprising  based  on  the  discussion 
from  Section  6.2  which  indicated  that  the  model-form  error  becomes  more  complex  for 
increasing  mode  numbers.  These  results  also  imply  that  the  PTeff  model  should  contain 
information  from  mode  6  (i.e.  the  next  asymmetric  mode)  as  well.  These  conclusions  indicate 
that  the  infonnation  from  a  sensitivity  analysis  can  be  used  to  guide  the  construction  of  the  PTejf 
model. 


Figure  6.17.  First-order  effect  sensitivity  indices  for  remaining  model-form  uncertainty  in 
effective  slope  model  as  a  function  of  Mach  number,  modes  1  to  4 


Co 


Figure  6.18.  Total  effect  sensitivity  indices  for  remaining  model-form  uncertainty  in 
effective  slope  model  as  a  function  of  Mach  number,  modes  1  to  4 
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6.5  Summary  and  Conclusions 

The  sources  of  model  uncertainty  were  investigated  for  aerodynamic  pressure  and  nonlinear 
structural  dynamic  response  predictions  in  the  coupled  aeroelastic  model  of  a  thin  panel 
subjected  to  hypersonic  flow.  The  model  uncertainty  was  quantified  through  an  internal  error 
model  at  the  local  slope  of  the  deformed  panel,  which  differs  from  other  approaches  that  add  an 
external  discrepancy  term  to  the  model  prediction.  The  benefit  of  capturing  the  model 
uncertainty  at  the  source  within  the  model  helps  maintain  the  underlying  physics,  while 
improving  the  accuracy  of  the  prediction.  This  approach  led  to  significant  error  reductions 
between  3ld-order  piston  theory  and  CFD  results,  as  seen  in  a  verification  study  performed  under 
steady  and  unsteady  flow  conditions  and  for  expected  panel  deformations.  The  solution 
approximation  error  in  the  structural  ROM  was  assessed  by  comparing  the  response  of  a  6-mode 
ROM  and  a  FEA  model,  both  coupled  with  piston  theory.  The  results  suggest  that  as  the  limit 
cycle  oscillation  (LCO)  amplitude  increases,  the  solution  approximation  errors  from  the  ROMs 
become  more  significant  in  the  structural  response.  The  aeroelastic  response  of  the  panel  was 
investigated  using  piston  theory  with  the  effective  slope-based  uncertainty  model  with  the  6- 
mode  ROM.  A  large  reduction  in  the  predicted  LCO  amplitude  was  observed  when  model-form 
uncertainty  was  accounted  for  in  the  aero-pressure  computations.  Finally,  sensitivity  analysis 
was  performed  to  assess  the  relative  significance,  on  the  LCO  amplitude,  of  the  structural  modes 
that  form  the  ROM  basis.  The  main  and  total  effect  sensitivity  indices  indicated  that  model 
uncertainty  in  the  pressure  component  from  mode  2  was  the  most  important  contributor  to  the 
total  variance.  Comparison  of  the  sensitivity  indices  for  modes  1,  3,  and  4  suggested  that 
interaction  plays  and  important  role  in  the  contribution  from  these  modes. 
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7  A  PRE-VALIDATION  STUDY  ON  LEGACY  HIGH-SPEED  WIND  TUNNEL  DATA 

One  of  the  greatest  challenges  to  the  development  of  scramjet-based  hypersonic  vehicles  is 
accurately  modeling  the  aerothermoelastic  response.  The  interaction  of  hypersonic 
aerodynamics,  structural  deformation,  and  heat  transfer  pose  multiple  analysis  challenges  [23]. 
The  computational  requirements  are  too  large  for  a  coupled,  direct  computational  fluid  dynamics 
and  finite  element  simulation  approach,  and  one  cannot  simply  factor-of-safety  out  of  the 
problem;  the  mass  constraints  to  meet  mission  requirements  are  too  tight.  Reduced-order  models 
are  needed,  and  the  errors  in  those  models  need  to  be  quantified  using  validation  data. 

Recent  studies  by  AFRL  and  its  affiliates  have  focused  on  a  reduced-order  model  approach 
to  panel  flutter  (one  of  the  structural  challenges  of  sustained  hypersonic  flight),  using  piston 
theory  as  the  basis  for  pressure  predictions.  This  group  has  even  pursued  validation  of  piston 
theory  using  published  wind  tunnel  data  originally  obtained  for  other  purposes  [9,15,23].  Their 
work  made  use  of  a  1986  study  conducted  by  Glass  and  Hunt  in  NASA’s  8’  High-Temperature 
Tunnel  (HTT).  While  the  study  by  Smarslok  and  Mahadevan  provides  a  framework  for  validation 
work,  it  lacks  an  essential  step  in  the  use  of  historical  data  for  validation  purposes. 

No  reduced-order  model  operates  on  its  own.  The  boundary  conditions  and  inputs  fed  into 
said  model  are  derived  from  predictions  produced  by  other  reduced-order  models.  For  example, 
in  validating  piston  theory  using  the  Glass  and  Hunt  1986  data,  one  must  establish  the 
“undisturbed”  post-shock  flow  on  which  the  instrumented  dome  impinges.  In  the  study  by 
Smarslok  and  Mahadevan  [23],  these  values  are  produced  using  2D  oblique  shock  theory.  If  2D 
oblique  shock  theory  is  not  appropriate,  or  if  there  are  biases  in  the  reported  data,  then 
inferences  about  the  validity  and  accuracy  of  piston  theory  drawn  from  the  Glass  and  Hunt  1986 
data  will  be  wrong.  It  is  to  avoid  this  type  of  error  that  a  pre-validation  study  is  conducted  and 
reported  here. 

7.1  Legacy  Aerothermal  Data  from  the  NASA  High-Temperature  Tunnel 

In  1986,  Glass  and  Hunt  studied  the  flow  over  a  shallow  dome  protuberance  on  a  flat  plate  in 
nominally  Mach  6.5  flow  to  investigate  the  thermal  and  structural  loads  on  body  panels  in  extreme 
environments  [44].  The  experiments  in  this  study  were  conducted  in  the  NASA  Langley  8’  HTT 
and  provide  rare  data  relevant  to  the  validation  of  piston  theory.  Originally,  these  data  were 
intended  to  shed  light  on  the  flow  over  shuttle  tiles  during  reentry. 


Figure  7.1.  Sketch  of  Apparatus  from  Glass  and  Hunt  1986  HTT  Experiments 

Domes  instrumented  with  pressure  taps  and  thermocouples  were  embedded  in  spherical  dome 
specimens  on  a  flat  plate  held  at  an  inward  five  degree  incline  with  respect  to  the  flow.  So, 
between  the  flow  and  the  dome  plate,  there  was  an  oblique  shock.  The  inviscid  flow  region  was 
still  supersonic  (M«5.7)  and  uniform  post-shock  until  encountering  the  spherical  dome 
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protuberances.  Figure  7.1  illustrates  the  structure  of  the  flow  studied.  There  were  pressure  data 
for  33  experiments  at  58  locations  on  the  dome. 

Given  the  simplicity  of  the  flow  and  the  ample  data  taken,  the  Glass  and  Hunt  1986  study  would 
seem  an  ideal  candidate  for  use  in  the  validation  of  piston  theory.  In  theory,  all  one  should  need  to 
do  is  compute  the  post-shock  flow  conditions  and  use  these  as  the  “undisturbed”  flow  conditions 
input  into  piston  theory.  The  deflections  induced  by  the  dome  are  simple  and  easy  to  compute. 
Therefore,  determining  the  discrepancy  between  piston  theory  prediction  and  real  pressures 
should  be  a  relatively  simple  task. 

However,  an  initial  review  of  the  data  reveals  that  the  ratio  of  the  measured  post-shock  static 
pressure  to  the  measured  freestream  pressure  does  not  match  the  static  pressure  ratio  obtained 
using  oblique  shock  theory.  The  Mach  numbers  reported  range  between  6.53  to  6.62.  The 
reported  deflection  angle  is  five  degrees.  These  conditions  should  produce  a  weak  oblique 
shocks  with  pressure  ratios  ranging  from  2.11  to  2.13.  However,  the  ratios  between  the 
measured  flat  plat  pressures  and  the  reported  freestream  pressures  range  from  1.85  to  2.00. 
Something  is  wrong:  a  bias,  either  in  the  data  or  in  its  theoretical  interpretation. 

The  pressure  discrepancy  is  important  in  that  it  throws  the  data,  the  interpretation  of  the  data,  or 
both  into  question.  If  the  post-shock  pressure  predicted  by  oblique  shock  theory  is  wrong,  then 
so  is  the  post-shock  Mach  number.  That  Mach  number  is  a  key  input  into  the  predictions  of  piston 
theory.  If  the  pressure  measurements  are  biased,  then  it  throws  the  pre-shock  and  post-shock  Mach 
numbers  into  question.  Finally,  if  the  discrepancy  is  due  to  3D  flow  effects,  those  same  effects 
will  impact  the  dome  pressure  measurements.  Whether  the  measurements  are  in  error  or  the 
theoretical  treatment  of  the  data  is  in  error,  there  is  an  error,  and  that  error  will  affect  any 
comparison  between  the  measurements  obtained  by  Glass  and  Hunt  and  the  predictions  of  piston 
theory.  The  identification  of  that  error  is  the  focus  of  this  investigation  and  constitutes  the  “pre- 
validation  study”  alluded  to  in  the  title. 

There  are  many  potential  physical  explanations  for  the  pressure  discrepancy.  Two  simple 
explanations  are  explored  here.  The  first  explanation  is  a  bias  in  the  reported  freestream 
pressures.  A  cursory  investigation  showed  that  a  10%  bias  in  freestream  pressure  could  readily 
account  for  the  discrepancy  in  pressure  ratios.  A  1973  study  of  the  facility  (NASA  8’  HTT) 
showed  that  the  pressure  profile  across  the  test  flow  does,  in  fact,  vary  on  the  order  of  10%  [46]. 
This  may  be  due  to  the  conical  nature  of  the  wind  tunnel  flow.  So,  if  the  freestream  pressure 
reported  was  drawn  from  a  different  part  of  the  flow  than  that  encountering  the  test  bed,  that 
would  certainly  explain  the  bias.  Alternatively,  if  the  deflection  (or  effective  deflection)  of  the  flow 
were  off  by  half  a  degree  from  the  nominal  five  degrees,  that  too  could  account  for  the  observed 
pressure  ratio  discrepancy.3 

Another  important  potential  explanation  for  the  pressure  ratio  discrepancy  is  pressure  relief  via 
a  strong  edge  vortex.  (Actually,  it  would  be  a  pair  of  edge  vortices,  one  on  each  side.)  This 
possibility  is  difficult  to  explore  quantitatively  because  the  authors  are  unaware  of  a 
computationally  efficient  edge  vortex  model.  In  fact,  even  CFD-based  approaches  to  supersonic 
vortex  modeling  are  known  to  suffer  from  substantial  numerical  dissipation.  Without  a  concrete 
model  structure  for  the  edge  vortex,  it  is  not  possible  to  statistically  test  this  potential  “bias” 
source. 

Finally,  pressure  measurements  on  the  dome  itself  are  highly  suggestive  of  an  edge  vortex.  In 
the  absence  of  an  edge  vortex,  pressure  taps  that  reflect  each  other  along  the  dome  diameter 
aligned  with  the  freestream  flow  should  have  equal  pressure  measurements.  To  be  sure,  there 
will  be  some  discrepancy  attributable  to  random  error,  but  at  first  glance,  the  scale  and 


71 

Approved  for  public  release;  distribution  unlimited. 


consistency  of  that  bias  would  seem  to  indicate  an  edge  vortex.  The  evidence  for  such  an  edge 
vortex  is  explored  on  that  basis  in  this  study. 

7.2  Mathematical  Methods 
7.2.1  Posterior  p-Value 

Assessing  the  hypothetical  sources  of  pressure  ratio  bias  presents  a  difficult  problem.  Given 
the  provisional  nature  of  these  hypotheses,  a  p-value  approach  seems  appropriate.  However,  these 
are  not  simple  hypotheses  corresponding  to  a  physical  model  with  a  set  of  parameters.  The  goal 
here  is  to  test  the  assertion  that  there  is  a  bias  in  the  freestream  pressure  measurements  or  in  the 
deflection  angle,  without  committing  to  the  magnitude  of  that  bias.  So,  it  is  a  test  of  model  form, 
rather  than  a  test  of  parameter  values. 

Readers  unfamiliar  with  traditional  p-values  may  benefit  from  a  quick  review.  In  short,  the 
“p”  in  “p-value”  stands  for  plausibility.  The  p-value  accorded  to  a  hypothesis  reflects  its 
plausibility  in  light  of  the  data.  Given  a  test  statistic,  the  p-value  for  a  hypothesis  is  the 
probability  of  having  gotten  a  less  favorable  test  score  given  that  the  hypothesis  is  true. 

That  is,  if  T  is  the  test-statistic,  and  small  values  of  T(x,  II)  are  favorable  to  hypothesis  H, 
then  the  p-value  for  H  is  as  follows:5 

Pls(#  I  x)  =  Prm  (T (x,H)  <  T(X',H)).  (7.1) 

If  large  values  of  T(x)  are  favorable  to  hypothesis  H,  then  the  p-value  of  H  is  as  follows: 

Pis (H  |  x)  =  ?rx)H  ( T(x,H )  >  T(X',H)).  (7.2) 

Given  a  test-statistic,  TJx),  the  user  detennines  what  kinds  of  values  are  considered  favorable  or 
unfavorable.  That  might  seem  arbitrary,  but  a  directional  preference  usually  suggests  itself.  For 
example,  suppose  /u  were  the  hypothetical  mean  of  X,  and 


T{x,v)  = 


k= 1 


(7.3) 


Here,  T  represents  a  discrepancy,  and  small  values  are  favorable  to  the  hypothesis  in  question. 
Alternatively,  suppose  instead  that  T  were  the  following  likelihood  ratio: 


T{x,0) 


f{x\0) 
max0,  f[x\6') 


(7.4) 


where /(x|  6)  is  the  probability  of  x  given  9.  In  this  second  example,  larger  values  of  T  are  more 
favorable  to  the  hypothesized  value  of  6.  While  not  all  test-statistics  are  as  simple  as  measuring 
goodness  or  badness  of  fit,  their  directional  interpretations  tend  to  remain  intuitive,  simply 
because  test  statistics  are  designed  to  support  a  simple  directional  interpretation. 

When  dealing  with  a  precise  hypothesis,  computing  a  p-value  is  relatively  straightforward.  If 
nothing  else,  one  can  generate  Monte  Carlo  replicates  under  the  assumption  that  said  hypothesis 
is  true  and  record  the  resulting  Monte  Carlo  sample  of  the  test  statistic.  The  p-value  is  the  number 
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of  test  statistic  replicates  of  the  test  statistic  that  are  worse  than  the  one  obtained  with  the  real 
data,  divided  by  the  total  number  of  Monte  Carlo  replicates. 

To  deal  with  the  compound  problem  of  simultaneously  tuning  the  model  parameters  and 
assessing  the  plausibility  of  the  model  form,  only  a  slight  variant  is  needed.  First,  tune  the 
parameters  to  yield  either  a  Bayesian  posterior  or  confidence  structure,  assuming  that  the 
hypothesized  model  form  is  correct.  Compute  the  test  statistic  for  the  observation  and  tuned 
parameter  values.  Next,  draw  a  Monte  Carlo  sample  of  the  parameters  from  the  posterior.  For 
each,  generate  a  random  value  set  of  the  data,  x;  re-infer  the  parameters  as  though  x  were  the 
original  observation;  and  recompute  the  test  statistic  value.  From  here,  the  calculation  of  the  p- 
value  is  the  same  as  for  a  precise  hypothesis.  However,  instead  of  testing  the  hypothesis 

x~f(x  \0)  (7.5) 

for  some  specified  8,  the  posterior  p-value  will  test  the  hypothesis  that 

3 e  s.t.  x~f(x\0 )  (7.6) 

Figure  7.2  demonstrates  the  power  of  this  approach  to  falsify  a  normality  hypothesis  (i.e.  that  a 
sample  has  a  normal  or  Gaussian  distribution)  when  the  data  are  in  fact  unifonn  or  Cauchy.  The 
test  statistic  used  is  explained  in  Section  7.2.3. 

P-values  and  posterior  p-values  provide  a  simple  plausibilistic  interpretation  for  statistical 
results.  A  high  p-value  reflects  that  the  hypothesis  in  question  is  still  plausible  (i.e.  not  yet 
falsified)  in  light  of  the  data,  and  a  low  p-value  reflects  the  hypothesis  in  question  is  now 
implausible  (i.e.  falsified)  in  light  of  the  data.  At  the  risk  of  redundancy,  it  should  be  stressed  that 
a  high  p-value  alone  does  not  prove  a  hypothesis  true.  In  fact,  given  a  set  of  data,  multiple 
competing  hypotheses  may  be  plausible.  It  is  only  when  one  has  eliminated  all  reasonable 
alternatives  that  a  given  hypothesis  may  be  taken  as  confirmed  by  the  data. 

Test  for  Normality,  Uniform  X  Test  for  Normality,  Cauchy  X 


Figure  7.2.  Falsification  Power  of  Posterior  p-Value  Approach  for  Various  Sample  Sizes 
(Light  Blue  =  10,  Dark  Blue  =  20,  Green  =  50,  Red  =  100) 
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7.2.2  Metropolis-Hastings  Algorithm 

Bayesian  inference  was  accomplished  in  this  study  using  a  slight  variant  on  the  Metropolis- 
Hastings  algorithm.  Metropolis-Hastings  does  not  perform  Bayesian  inference  exactly;  it  yields 
a  Monte  Carlo  sample  distributed  according  to  a  specified  likelihood  function:  in  this  case,  the 
posterior.  The  variant  used  here  involves  working  with  iteration  on  a  large  sample,  rather  than 
selection  from  a  recorded  Markov  Chain.  What  this  variant  lacks  in  computational  efficiency,  it 
makes  up  for  in  simplicity  and  robustness. 

The  process  is  simple.  Start  with  a  seed  sample,  x0,  of  desired  size.  A  size  of  10,000  is 
sufficient  for  most  Monte  Carlo  analyses.  Next,  iterate  as  follows: 

1 .  Add  a  random  perturbation,  ek,  to  current  sample,  x*.  Record  the  result,  xl  =  x/f  + 
tk-  The  perturbation,  e*,  must  be  an  i.i.d.  sample  of  a  distribution  that  is  symmetric 
about  zero. 


2. 

3. 

4. 


Take  the  ratio  of  likelihoods,  R  =  — 7 — - 

Take  a  random  selection  seed,  u,  uniformly  distributed  on  [0,  1], 

Assign  members  of  the  current  sample  and  the  perturbed  sample  to  the  new  sample 
according  to  whether  the  ratio  of  a  member  of  the  sample  is  greater  than  the 

selection  seed.  That  is,  xk+li  =  [[f?,  >  +[[f?;  <ui]jxk  i,  where  [[  ]]  returns  a 


Boolean  truth  value  of  zero  or  one.  This  formula  results  in  a  probability  of  //,  or  1 
(whichever  is  less)  of  replacing  Xkj  with  x  . 

Continue  this  process  until  the  empirical  CDF  of  the  Monte  Carlo  sample  has  stopped  making 
significant  changes. 

When  using  the  algorithm  outlined  above,  the  user  should  be  mindful  of  the  following  details: 

•  The  perturbation  distribution  should  be  symmetric  or  near-symmetric.  A  normal 
distribution  with  zero  mean  is  a  safe  option,  so  long  as  non-physical  samples  (if  any 
occur)  are  rejected. 

•  It  helps  to  have  the  perturbation  distribution  be  on  the  scale  of  the  target  distribution. 
The  user  may  or  may  not  be  able  to  guess  what  that  scale  is. 

•  Since  the  likelihood  ratio  drives  the  selection  algorithm,  the  user  does  not  need  to 
know  the  normalization  coefficient  of  the  distribution  you  are  sampling.  This  is  one 
of  the  most  useful  features  of  Metropolis-Hastings. 

•  Another  major  advantage  of  Metropolis-Hastings  over  other  methods  is  its 
applicability  to  multidimensional  problems.  There  is  nothing,  in  principle,  to 
distinguish  a  likelihood  over  one  variable  from  a  likelihood  over  several. 

•  However,  the  larger  the  inference  problem  (in  terms  of  dimensionality),  the  more 
iterations  it  takes  to  converge  Metropolis-Hastings. 

A  short  Matlab  example  may  be  instructive.  Suppose  the  analyst  wants  to  sample  from  a 
gamma  distribution  with  scale  parameter  =  1  and  shape  parameter  =  1/2.  The  probability 
density  function  for  this  variable  is  as  follows: 


f(x)  =  x  2e~x 


(7.7) 
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The  analyst  could  map  from  a  uniform  random  sample  to  this  distribution  by  mapping  with  the 
inverse  of  the  incomplete  gamma  function.  Or,  alternatively,  the  analyst  could  generate  a  sample 
more  easily  using  the  Metropolis-Hastings  algorithm. 

Let  N  be  the  desired  size  of  the  sample.  The  following  Matlab  script  will  generate  an  i.i.d. 
sample  of  size  N  from  the  gamma  distribution  with  unity  scale  and  shape  of  one  half: 

xO  =  -  sqrt(2)  *  log  (  1  -  rand(N,l)  ); 

%  initial  guess  from  exponential  distr 

x  =  xO; 

for  k  =  1:100 

x  pr  =  x  +  sqrt(2)  *  randn (  N  ,  1  );  %  normal  perturbation 
R  =  exp (  x  -  x  pr  )  . *  sqrt (  x  . /  x  pr  )  . *  (  x  pr  >=  0  ) ; 
u  slot  =  rand (  N  ,  1  ) ;  %  uniform  sample  for  selection  x  =  x 
pr  . *  (  R  >=  u  slct  )  +  x  .*  (  R  <  u  slct  ); 

%  final  selection  of  new  xs 
end 

How  U  is  found  U  Pool  Defect  Area 


Figure  7.3.  Calculation  of  the  U-Pool  Defect  Statistic 
7.2.3  U-pool  Defect 

The  U-pool  defect  is  a  simple  and  widely  applicable  test  statistic.  It  was  originally  suggested 
by  Ferson  and  Oberkampf  in  their  response  to  the  Sandia  Challenge  Problem  and  makes  use  of  the 
uniformity  of  CDF  values  [70].  Suppose  one  has  a  sample,  x,  and  it  is  assumed  that  there  is  some 
value  of  9  such  that  x  ~  fix\9).  Then,  for  some  future  sample,  x,  the  CDF  of  the  predictive 
distribution  based  on  x  should  be  uniformly  distributed.  That  is  F(x]x)  should  be  an  i.i.d. 
unifonn  sample.  The  U-pool  defect  is  the  area  of  the  absolute  difference  between  the  empirical 
CDF  of  F{x  |x)  and  the  line  x  =y  over  [0,1]. 

In  this  study,  a  slight  variant  of  the  U-pool  defect  is  used;  this  variant  could  be  called  the  Auto- 
U-pool  defect.  Instead  of  taking  a  second  sample,  xf,  x  is  simply  fed  back  in  as  though  it  were  a 
new  sample.  So,  strictly  speaking,  F(x]x)  with  xl  =  x,  won’t  be  i.i.d.  uniform,  but  as  long  as  the 
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sample  size  is  sufficiently  larger  than  the  number  of  parameters,  it  will  be  close  enough.  Moreover, 
for  small  sample  sizes,  Auto-U-pool  defect  values  tend  to  be  smaller  than  they  would  be  for  an 
actual  i.i.d.  uniform  sample;  so,  it  will  tend  to  give  high  p-values,  i.e.  conservative  inference. 

7.2.4  Confidence  Structure  on  the  Non-Parametric  Difference 

Confidence  structures  are  a  simple  and  robust  alternative  to  Bayesian  inference  [71].  There  are 
three  differences  between  a  confidence  structure  and  a  Bayesian  posterior.  First,  a  confidence 
structure  does  not  require  a  prior  distribution.  Second,  while  a  Bayesian  posterior’s  belief  values 
usually  have  a  subjective  interpretation,  a  confidence  structure’s  belief  values  have  a  statistical 
confidence  (i.e.  coverage  probability)  interpretation.  Finally,  while  Bayesian  belief  values  satisfy 
the  Kolmogorov  axioms,  confidence  structures  only  satisfy  the  more  general  Shafer  axioms. 
This  generality  enables  the  application  of  confidence  structures  to  problems  of  non-parametric 
statistics,  which  is  useful  in  this  study. 

The  Mann- Whitney  U-statistic  is  the  foundation  of  a  well-known  test  for  the  equality  of  two 
distributions  given  no  assumptions  about  the  form  of  those  distributions.  Moreover,  under  the 
assumption  that  the  two  distributions  considered  are  identical  except  for  an  offset,  Mann-Whitney 
U  is  a  pivot.  That  is,  U(x  -  8,y ),  where  8  is  the  offset  between  the  distributions  of  x  and  y,  has  a 
known  distribution.  The  pivot  is  the  key!  It  has  previously  been  demonstrated  that  a  pivot  can  be 
used  to  readily  construct  a  confidence  structure  [71].  The  confidence  structure  developed  from 
the  Mann-Whitney  U-statistic  is  referred  to  in  this  work  as  the  “Non-Parametric  Difference.” 

Figure  7.4  illustrates  how  a  pair  of  samples  translates  to  a  non-parametric  difference.  First, 
one  takes  all  m  x  ni  possible  values  ofxi  —xi  and  sorts  them  in  ascending  order.  This  set  of  points 
partition  of  the  real  line.  This  partition  is  the  set  of  focal  elements  for  the  non-parametric 
difference.  Numbering  these  focal  elements,  Ak,  from  0  to  m  x  n2,  the  weight  assigned  to  Ak  is 
equal  to  the  probability  that  U(x  -  8,y )  =  k,  where  8  is  the  true  unknown  offset. 


Empirical  CDFs  for  Two  Samples  Non-Parametric  Difference  Derived  from  Above  Samples 


Figure  7.4.  Example  of  Non-Parametric  Difference  Confidence  Structure 


7.3  Data  and  Models 

Three  of  the  four  potential  explanations  for  the  freestream  bias  are  explored  using  the  posterior 
p-value  method  described  in  Section  7.2.1.  Each  of  these  hypotheses  leaves  several  free 
parameters.  First,  it  is  assumed  that,  even  in  the  absence  of  bias,  there  is  still  a  random  error  in 
the  measurements  of  the  pre-shock  and  post-shock  pressures.  The  random  error  is  assumed  to 
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be  normal  (i.e.  Gaussian).  However,  the  standard  deviation  of  these  normal  errors  are  not 
known.  So,  in  addition  to  whatever  bias  tenn  is  inferred,  the  true  pre-shock  pressure  in  each  flow 
is  inferred,  as  are  the  standard  deviations  of  the  errors  in  the  pre-shock  and  post-shock  pressure 
measurements.  The  following  subsection  describe  the  exact  mathematical  models  used  to 
support  these  inferences. 

7.3.1  Pre-shock  Conditions 

For  the  analyses  in  Sections  7.2,  it  is  assumed  that  there  is  some  error  (bias,  random,  or  both)  in 
the  reported  pressures.  As  mentioned  above,  the  true  unknown  pre-shock  pressure  is  inferred  via 
Bayesian  inference.  To  the  authors’  knowledge,  there  is  no  such  thing  as  a  pure  Mach  number 
measurement  device.  Moreover,  there  is  no  record  of  direct  velocity  measurement  (e.g.  via  laser 
doppler  anemometry)  in  the  Glass  and  Hunt  1986  report.  So,  it  is  reasonable  to  assume  that  the 
Mach  numbers  reported  in  Glass  and  Hunt  [44]  were  based  on  the  pre-shock  pressure 
measurements,  as  described  in  the  1973  wind  tunnel  calibration  report  [46].  However,  if  the 
reported  pressure  has  an  error,  so  does  the  reported  Mach  number. 

Inferred  values  of  the  true  pre-shock  pressure  imply  tme  pre-shock  Mach  numbers  that  are 
(slightly)  different  from  the  reported  Mach  number.  It  is  assumed  that  the  reported  pre-shock 
Mach  number,  M»,  corresponds  correctly  to  the  reported  pre-shock  pressure  value,  px.  The 
corrected  Machnumber  is  adjusted  to  the  inferred  pressure  on  the  assumption  that  the  entropy  is 
near  equal  between  reported  and  corrected  conditions.  This  leads  to  the  following 
transformation: 


f 

y— 1  ^ 

2  ^ 

f  y  —  l  ,3 

i  i  '  \/i 2 

(  P  ^ 

7 

-1 

l  2  "J 

\Pkt  ) 

V  J 


(7.8) 


where  y  is  the  ratio  of  specific  heats  (approximately  equal  to  1.4  at  test  section  conditions), 
pJf  is  the  inferred  true  pre-shock  pressure  value,  and  is  the  resulting  corrected  pre-shock 
Mach  number.  Using  this  information,  the  true  value  of  the  post-shock  (i.e.  flat  plate)  pressure  can 
be  calculated  as  well,  according  the  bias  models  described  in  the  next  three  subsections. 

In  the  bias  models  explored,  it  is  assumed  that  the  2D  oblique  shock  relations  apply  and  that 
either  random  error,  bias  in  the  freestream  measurements,  or  bias  in  the  deflection  were 
responsible.  Under  these  assumptions,  the  relationship  between  the  true  pre-shock  (p,T )  and 
true  post-shock  (pfPT)  pressures  is  given  as  follows: 


lML  =  \  +  2^- 


p^t  r+ 

where  is  the  shock  angle,  calculated  as  follows: 
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where  6  is  the  deflection  angle  (i.e.  the  angle  by  which  the  flow  is  turned  inward)  and  X  and  %  are 
parameters  defined  as  follows: 


(7.11) 


(7.12) 


This  explicit  fonnulation  for  the  oblique  shock  relations  is  available  in  a  commonly  used  text  on 
compressible  aerodynamics  [72]. 

7.3.2  No-Bias  Hypothesis 

Under  the  “no  bias”  hypothesis,  it  is  assumed  that  there  is  no  bias  in  the  measurement  of 
either  the  pre-shock  or  post-shock  pressure,  i.e.  that  any  random  is  error.  That  is,  the  observed 
bias  between  the  predicted  and  observed  pressure  ratios  across  the  shock  is  explained  as  a 
random  fluke. 

Given  sample  values  for  the  true  values  of  pre-shock  and  post-shock  pressures,  it  is  possible 
to  calculate  the  posterior  likelihood  of  said  sample.  This  is  a  necessary  step  in  executing 
Bayesian  inference  to  tune  the  parameters  in  the  different  bias  models.  In  each,  a  flat  prior  for  the 
true  pre-shock  pressures  is  used.  Moreover,  a  i  prior  is  used  for  the  two  unknown  standard 

deviations.  Given  the  normality  assumptions  described  above,  the  likelihood  can  be  computed  as 
follows: 
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where  k  indexes  the  runs,  OfP  is  the  standard  deviation  of  the  random  error  in  the  post-shock  (flat 
plate)  pressure  measurements,  and  is  the  standard  deviation  of  the  random  error  in  the  pre¬ 
shock  pressure  measurements.  Both  of  these  standard  deviations  are  unknown;  so,  n  +  2 
parameters  are  being  inferred.  There  are  2 n  pressure  observations  (one  on  each  side  of  the  shock 
for  each  run);  so,  as  long  as  more  than  two  runs  are  used,  the  problem  is  closed  (i.e.  detenninate).  It 
is  this  closure  that  allows  the  use  of  non-infonnative  priors. 

7.3.3  Free  Stream  Bias 

Under  the  “freestream  bias”  hypothesis,  it  is  assumed  that  there  is  a  proportional  bias  in  the 
measurements  of  freestream  pressure,  in  addition  to  random  error.  The  resulting  likelihood  model 
is  as  follows: 
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where  AT  is  the  proportional  bias,  which  is  an  additional  inferred  parameter.  Interestingly,  this 
bias  is  equivalent  to  having  a  proportional  bias  in  the  measured  total  (i.e.  chamber)  pressure, 
meaning  that  this  model  covers  two  (indistinguishable,  given  the  data)  possible  bias  explanations. 

7.3.4  Deflection  Bias  Hypothesis 

Under  the  “deflection  bias”  hypothesis,  it  is  assumed  that  the  deflection  is  other  than  the 
reported  nominal  5  degrees.  The  pressure  values  and  likelihoods  are  related  as  expressed  as  in 
Equations  1-2,  the  deflection  angle  is  merely  different.  This  difference,  §6,  is  the  tunable  bias 
parameter. 

7.3.5  Dome  Pressure  Asymmetry 

If  the  pressure  bias  is  due  to  relief  from  an  edge  vortex,  this  edge  vortex  will  also  cause 
pressure  to  decrease  from  the  center  towards  edges.  Three  cases  were  examined:  7in,  14in, 
28in.  The  relative  position  of  the  pressure  dome  in  these  three  cases  is  illustrated  in  Figure 
7.5.  Fortunately,  in  the  7in  and  14in  runs,  the  pressure  dome  is  located  off-center.  Moreover,  on 
the  front  half  of  the  pressure  dome,  each  tap  has  another  tap  reflected  across  the  centerline 
oriented  with  the  freestream  flow  direction.  A  span-wise  decrease  in  pressure  caused  by  an 
edge  vortex  would  show  up  as  an  asymmetry  in  the  pressures  between  paired  pressure  taps. 
Conversely,  in  the  absence  of  an  edge  vortex,  the  paired  pressure  taps  should  produce 
measurements  that  are,  on  average,  equal. 

The  asymmetry  of  paired  dome  pressures  was  explored  using  the  non-parametric  difference 
described  in  Section  7.2.4.  It  is  assumed  that  each  pair  of  taps  has  the  same  random  error 
distribution.  It  is  further  assumed  that  the  underlying  offset  between  a  given  pair  of  pressure  taps 
has  a  constant  value  throughout  all  runs,  obscured  only  by  random  experimental  error.  The  non- 
parametric  difference  may  be  taken  as  a  direct  measure  of  the  asymmetry  in  each  pair. 

P-values  derived  from  the  Mann-Whitney  U-test  also  provide  a  holistic  assessment  of 
asymmetry  across  pressure  dome.  These  p-values  were  split  according  to  dome  diameter  (7in, 
14in,  or  28in).  If  the  underlying  dome  pressures  were  symmetric,  then  the  Mann-Whitney  p- 
values  should  be  uniformly  distributed.  The  U-pool  defect  of  Section  7.2.3  was  then  used  to 
assess  whether  or  not  they  were,  and  whether  their  deviation  from  uniformity  was  statistically 
significant. 

7.4  Results 

7.4.1  Pressure  Ratio  Bias 

The  methods  applied  in  this  study  were  chosen  for  their  robustness,  not  their  computational 
efficiency.  As  such,  it  was  only  possible  to  generate  100  posterior  bootstrap  samples.  Moreover, 
the  limited  samples  still  provided  estimates  of  the  p-values  for  the  different  hypotheses.  Those 
results  are  as  follow: 

-  No  bias  hypothesis:  p-value  =  0  out  of  1 00  (tme  pis  <  3%  with  95%  confidence) 

Freestream  bias  hypothesis:  p-value  =  24  out  of  100  (true  pis  between  16%  and  34%  with 
95%  confidence) 
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Deflection  bias  hypothesis:  p-value  =  15  out  of  100  (true  pis  between  9%  and  24%  with 
95%  confidence) 

None  of  the  three  hypotheses  examined  using  the  posterior  p-value  approach  merited  a  high 
plausibility.  However,  only  the  “no  bias”  hypothesis  was  conclusively  falsified.  It  bears  noting  that 
the  mismatch  statistic  (i.e.  the  U-pool  defect)  value  obtained  for  the  “bias”  case  is  so  much  larger 
than  its  bootstrap  counter-parts  that  using  any  reasonable  curve-fit  to  those  bootstrap  samples 
yields  a  machine-zero  p-value  for  the  “no  bias”  hypothesis.  That  is  to  say,  the  p-value  of  the  no¬ 
bias  hypothesis,  when  approximated  using  a  curve  fit  to  the  Monte  Carlo  statistic  sample  values, 
is  smaller  than  the  computer  is  able  to  represent. 

7.4.2  Pressure  Tap  Asymmetry 

Using  the  aggregated  Mann- Whitney  p-values,  it  was  found  that: 

In  the  7in  case,  symmetry  is  43.5%  plausible; 

In  the  Min  case,  symmetry  is  0.67%  plausible; 

In  the  28in  case,  symmetry  is  30.7%  plausible. 

Symmetry  is  statistically  implausible  for  the  14  inch  case,  where  the  pressure  dome  is  next 
to  the  edge,  statistically  plausible  for  the  28  inch  dome  case,  in  which  the  dome  is  centered,  and 
plausible  for  the  7  inch  dome  case  in  which  the  dome  is  not  centered,  but  also  not  close  to  the 
edge.  These  results  are  consistent  with  the  existence  of  a  significant  edge  vortex. 

Figure  7.6  illustrates  the  level  of  asymmetry  at  three  point-pairs  in  the  14-inch  case.  The 
asymmetry  is  stronger  for  points  more  widely  separated.  This  too  is  consistent  with  the 
existence  of  a  significant  edge  vortex.  The  front-most  point-pair  has  a  backwards  (negative) 
offset  relative  to  the  other  two  pairs.  It  is  unknown  whether  the  effect  is  due  to  a  real  underlying 
physical  phenomenon  (e.g.  shock-vortex  interaction)  or  to  random  measurement  error. 

28  inch  dome  14  inch  dome  7  inch  dome 


Figure  7.5.  Diagram  of  spherical  dome  specimens  and  positions 
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Figure  7.6.  Non-Parametric  Difference  Across  Three  Pressure  Tap  Pairs 
7.5  Conclusions  and  Future  Work 

This  research  presented  a  pre-validation  study  of  legacy  aerodynamic  data  obtained  in  1986  by 
Glass  and  Hunt  in  the  NASA  8’  HTT.  An  apparent  bias  in  the  pressure  data  confounds  the  direct 
application  of  these  data  to  the  validation  of  piston  theory.  Advanced  statistical  methods  were 
applied  to  explore  this  bias. 

It  is  unclear  how  the  observed  pressure  bias  in  the  Glass  and  Hunt  data  originates,  but  it  is 
clear  that  the  bias  is  the  result  of  unaccounted  physics  or  apparatus  bias,  and  not  an  artifact  of 
random  coincidence.  Furthermore,  there  is  positive  evidence  for  an  edge  vortex  in  the  Glass  and 
Hunt  data.  Whether  that  edge  vortex  is  solely  responsible  for  the  discrepancy  between  expected 
and  observed  ratios  of  pre-shock  and  post-shock  pressures  remains  unknown.  Moreover,  there  may 
be  other  complex  phenomena  that  would  affect  the  comparison  of  reduced-order  model 
predictions  to  the  Glass  and  Hunt  data. 

While  the  conclusions  presented  here  are  limited,  the  methods  used  in  this  study  should  be 
new  to  many  readers.  The  posterior  p-value  enables  the  assessment  of  hypotheses  about  model 
form,  even  when  model  parameters  are  left  free.  The  Metropolis-Hastings  algorithm  provides  a 
robust  way  to  accomplish  the  model  tuning  involved  in  both  traditional  Bayesian  inference  and  the 
more  computationally  intense  posterior  p-value  method.  In  a  different  vein,  the  non-parametric 
difference  allows  a  specialized  analysis  free  from  assumptions  about  distribution  form.  These 
are  the  kind  of  detailed  statistical  tools  that  engineers  will  need  if  they  are  to  undertake  verification, 
validation,  and  uncertainty  quantification  successfully. 

Pre-validation  work  with  the  Glass  and  Hunt  data  is  not  complete.  It  is  possible  that  more 
severe  statistical  tests  (e.g.  Kolmogorov-Smimov)  may  provide  a  basis  on  which  to  falsify 
hypotheses  explored  in  this  study.  Perhaps  more  important  is  the  problem  of  the  edge  vortex.  It 
is  impossible  to  account  for  the  effects  of  an  edge  vortex  without  a  coherent  hypersonic  edge 
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vortex  model.  To  the  authors’  knowledge,  this  remains  one  of  the  outstanding  problems  in 
aerodynamics. 

While  there  will  certainly  be  some  point  at  which  the  analyst  must  say  “enough  is  enough,” 
right  now,  the  unexplained  discrepancies  in  the  Glass  and  Hunt  data  are  too  large  to  ignore. 
Moreover,  it  is  anticipated  that  legacy  hypersonic  wind  tunnel  data  are  replete  with  such 
discrepancies.  The  wind  tunnel  engineers  of  previous  decades  simply  did  not  design  their 
experiments  with  quantitative  model  validation  in  mind.  However,  judicious  pre-validation 
studies,  like  the  one  commenced  here,  may  ultimately  provide  a  framework  for  using  those  legacy 
data  in  today’s  validation  efforts. 
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8  UNCERTAINTY  QUANTIFICATION  OF  STATE  BOUNDARIES  IN  THIN  BEAM 
BUCKLING  EXPERIMENTS 

State  selection  -  or  model  selection  -  has  garnered  interest  in  many  fields  including 
phylogenetics  [73],  social  research  [74],  ecology  and  evolution  [75],  protein  evolution  [76],  and 
data  mining  [77].  Selection  of  the  best  model  is  based  either  on  classical  hypothesis  testing, 
infonnation  theoretic  measures,  or  Bayesian  hypothesis  testing.  Classical  approaches  test 
likelihood  ratio  statistics  to  find  the  best  model  [78].  Infonnation  theoretic  measures  are  the 
Akaike  information  criterion  (AIC)  [79]  and  its  generalization  Network  information  criterion 
(NIC)  [80],  conected  Akaike  information  criterion  (AICc)  [81],  the  focused  infonnation  criterion 
(FIC)  [82],  the  Kullback-Leibler  information  criterion  (KLIC)  [83],  Takeuchi  information 
criterion  (TIC)  [84],  minimum  description  length  (MDL)  [85],  and  deviance  infonnation 
criterion  (DIC)  [86].  The  popular  Schwarz  Bayesian  infonnation  criterion  (BIC)  [87]  is  not  an 
infonnation  theoretic  measure  and  was  originally  derived  to  select  the  posterior  most  probable 
model.  BIC  is  often  taken  as  an  approximation  to  Bayes  factors  [88],  which  are  an  element  of  the 
Bayesian  posterior  odds  ratio.  Bayes  factors,  and  by  extension,  posterior  odds,  are  a  natural 
implementation  of  Occam’s  razor  [89],  and  require  no  modifications  for  model  complexity.  No 
assumptions  about  nested  models  or  estimates  of  model  complexity  are  required,  as  with  many 
non-Bayesian  methods.  Further,  they  are  suitable  for  comparing  multiple  hypotheses  (i.e. 
models),  where  most  classical  statistical  methodologies  only  compare  two  hypotheses. 
Computation  of  a  Bayesian  posterior  odds  [90]  does  require  prior  distributions  for  parameters  (or 
prior  probabilities  for  models)  to  be  assumed,  which  can  be  a  difficult  task.  Quantifying  model 
selection  uncertainty  is  straightforward  from  the  Bayesian  posterior  odds  ratio,  which  yields 
posterior  probabilities  for  each  model  which  can  then  be  used  to  rank  the  models[88,91]. 
Bayesian  model  averaging  [92],  which  mitigates  overconfidence  in  the  model,  follows  naturally. 
Note  that  estimation  of  model  selection  uncertainty,  however,  is  also  possible  using  model 
selection  criterion  such  as  AIC  or  BIC.  These  criteria  compute  model  weights  that  can  be 
interpreted  as  probabilities. 

State  boundary  uncertainty  quantification  (SBUQ)  quantifies  the  uncertainty  in  a  boundary 
between  two  states,  as  opposed  to  the  uncertainty  in  model  selection.  SBUQ,  of  course,  requires 
knowledge  of  the  model  selection  uncertainty.  Previous  work  on  this  problem  includes  that  of 
Hombal  and  Mahadevan  [98],  who  developed  a  probabilistic  treatment  of  model  selection  using 
error  surrogates  that  considers  uncertainty  in  the  classification  boundary.  Alternative 
fonnulations  of  the  SBUQ  problem  may  be  possible,  each  with  their  own  pros  and  cons. 
Stochastic  optimization  can  provide  a  distribution  of  the  critical  boundary  value  while  being  a 
simpler  overall  methodology,  but  do  not  provide  a  direct  means  for  belief  updating  in  the 
presence  of  new  information.  Discriminative  classification  techniques  (e.g.,  logistic  regression 
[99],  neural  networks  [100],  support  vector  machines  [101],  relevance  vector  machines  [102]) 
may  provide  a  means  of  obtaining  an  estimate  of  the  discrete  state  with  an  associated  confidence 
by  regressing  on  discrete  and  continuous  system  inputs.  Generative  classification  techniques 
(e.g.,  Gaussian  discriminative  analysis  [103])  compare  likelihoods  between  competing  models  to 
estimate  the  discrete  state.  However,  individually,  either  style  of  classification  may  not  account 
for  enough  information  to  confidently  determine  the  state.  Discriminative  methods  may  assume 
an  incorrect  causal  interpretation  of  the  system  (i.e.,  responses  do  not  cause  the  discrete  state) 
and  pose  difficulties  when  the  regressor  variables  are  not  independent.  Generative  methods 
require  stronger  assumptions  about  the  underlying  distribution  of  the  data  than  discriminative 
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methods  and  cannot  account  for  uncertainty  in  the  training  data  labels,  making  predictions 
overconfident. 

The  SBUQ  problem  was  addressed  by  expressing  the  Bayesian  posterior  odds  model  selection 
formulation  of  Carlin  and  Chib  [93]  and  Green  [94]  is  expressed  using  a  Bayesian  network  (BN) 
[95].  The  approach  naturally  provides  posterior  probabilities  for  each  candidate  model,  thus 
providing  information  about  the  uncertainty  in  model  selection.  The  models  under  comparison 
are  generally  non-nested  and  the  number  of  parameters  may  be  difficult  to  quantify,  particularly 
if  surrogate  models  or  hierarchical  Bayesian  models  are  used  to  better  quantify  model 
uncertainty.  The  BN  represents  a  joint  probability  distribution  in  factored  form.  BNs  are  well 
qualified  for  uncertainty  quantification  tasks  due  to  their  probabilistic  nature.  Their  graphical 
representation  facilitates  system  understanding  and  makes  independence  assumptions  clear.  The 
trained  BN  can  then  be  invoked  to  perform  model  selection  uncertainty  quantification  and 
estimate  the  uncertainty  of  the  critical  boundary  values.  At  the  cost  of  being  more  complex  than 
the  alternatives,  the  methodology  developed  in  this  study  combines  discriminative  and  generative 
classification  methods  within  the  Bayesian  network,  providing  a  logical  system  description  that 
accounts  for  uncertainty  in  the  state  boundary  and  allows  for  belief  updating  given  new 
information. 

Section  8.1  describes  a  physical  example  of  SBUQ  in  quantifying  the  uncertainty  in  the 
buckling  temperature  of  a  thin  beam.  Section  8.2  explains  the  SBUQ  methodology,  including 
construction  of  the  Bayesian  network,  handling  uncertainties  in  training  data,  and  quantifying  the 
uncertainty  in  the  state  boundary.  Section  8.3  uses  data  collected  from  thin  beam  buckling 
experiments  to  demonstrate  SBUQ  by  quantifying  the  uncertainty  in  the  buckling  temperature  of 
a  thin  beam,  which  represents  the  boundary  between  the  pre-buckled  and  post-buckled  states. 

8,1  Physical  Example:  Thin  Beam  Buckling  Temperature 

Shukla  and  Mignolet  [104]  performed  experiments  on  11  flat  beam  specimens  for 
identification  of  the  uncertainty  in  nonlinear  reduced  order  models  (ROMs).  Variance  in  the 
experimentally  detennined  natural  frequencies  was  attributed  to  the  uncertainty  in  the  preload 
induced  by  the  boundary  conditions  and  material  parameters.  However,  Perez  et  al.  [36],  while 
using  the  data  from  Shukla  and  Mignolet  [104]  to  calibrate  a  finite  element  analysis  (FEA),  noted 
the  possibility  that  the  flat  beam  specimens  had  buckled.  The  original  experiments  were  not 
looking  for  buckling,  but  the  slight  changes  that  buckling  causes  to  the  shape  of  the  slender  beam 
could  have  easily  gone  unobserved.  Unfortunately,  the  amount  of  experimental  data  was 
insufficient  to  prove  buckling  was  occurring  and  partially  responsible  for  the  variance  of  the 
experimental  data.  Figure  8.1  illustrates  the  inconclusiveness  of  the  buckling  hypothesis.  The  test 
data  point  (red  ‘+’)  in  the  lower  left  of  Figure  8.1  could  be  part  of  either  1)  a  ROM  that  does  not 
consider  buckling  (red  line,  quadratic  fit  of  data)  or  2)  FEA  (blue  lines)  that  considers  buckling. 
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Figure  8.1.  Shukla  and  Mignolet  [104]  experimental  results  compared  to  finite  element 
analysis  for  pre-  and  post-buckled  operating  states. 

In  an  attempt  to  determine  which  explanation  was  best  and  create  a  better  model,  a  new 
experiment  for  the  flat  beam  was  designed.  The  focus  of  these  tests  was  to  reduce  uncertainty 
induced  by  the  boundary  conditions  and  provide  sufficient  data  for  calibrating  the  FEA  and 
determining  the  presence  of  buckling.  This  was  achieved  by  using  a  single  beam  specimen  that 
remained  clamped  in  place  throughout  the  experiment.  The  beam  was  repeatedly  heated  until  a 
steady  state  temperature  was  reached,  then  it  was  allowed  to  cool.  Measurements  were  obtained 
as  the  beam  cooled.  While  distinct  pre-  and  post-buckled  behavior  was  not  the  focus  of  the 
Shukla  and  Mignolet  [104]  experiments,  Perez  et  al.  [36]  investigated  and  confirmed  this 
phenomenon.  The  case  study  presented  in  Section  8.3  uses  this  data  to  quantify  the  uncertainty  in 
the  buckling  temperature  using  the  state  boundary  uncertainty  quantification  methodology. 

The  theoretical  Euler  buckling  load  of  a  column,  Fcr,  may  be  calculated  using  classical 
mechanics  by 


F  — 

1  cr 


(8.1) 


where  E  is  the  modulus  of  elasticity,  /  is  the  minimum  moment  of  inertia  of  the  cross  sectional 
area,  and  Le  is  the  effective  length  of  the  column  considering  the  boundary  conditions.  For  a 
column  with  both  ends  fixed,  the  effective  length  is  Le  —  0.5 L,  where  L  is  the  actual  column 
length.  For  this  research,  critical  buckling  temperature  is  of  interest.  The  force  induced  by  a 
temperature  load  on  a  beam  fixed  at  both  ends  is 


Ft  =  —EAa(T0  -  T) 


(8.2) 
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where  a  is  the  coefficient  of  thermal  expansion,  T0  is  the  reference  temperature  where  the  beam 
has  zero  stress  and  zero  strain,  and  T  is  the  current  temperature.  From  Eq.  (8.1)  and  Eq.  (8.2)  and 
with  Le  —  0.5 L,  the  buckling  temperature  of  a  fixed-fixed  beam  is  shown  in  Eq.  (8.3). 

4tt2/  3, 

Ttuct  =  +  To  <S3> 

Alternatively,  buckling  may  be  considered  to  occur  when  the  tangent  stiffness  of  the  beam 
tends  towards  zero.  At  this  point,  the  beam  buckles  because  it  cannot  tolerate  any  additional  axial 
load.  As  a  result,  one  or  more  of  the  eigenvalues  have  a  zero  value,  and  the  natural  frequencies 
corresponding  to  those  eigenvalue  are  also  zero.  The  result  is  bifurcation  in  the  natural 
frequencies  at  the  buckling  temperature.  A  detailed  discussion  is  available  from  Kosmatka  [105]. 
Figure  8.2  shows  a  plot  of  the  first  and  third  natural  frequencies  of  the  thin  beam  as  computed  by 
an  idealized  finite  element  model.  The  bifurcation  in  the  first  and  third  natural  frequencies  occurs 
at  the  buckling  load,  x-axis  is  AT  —  T0  —  T. 


AT  AT 

Figure  8.2.  Idealized  Buckling  Behavior 

An  idealized  finite  element  analysis  of  the  thin  beam  produces  these  plots  of  first  (left)  and 
third  (right)  natural  frequencies  vs  temperature. 

The  idealized  behavior  indicates  that  a  precise  buckling  temperature  is  identifiable.  However, 
in  practice  this  is  not  the  case,  as  many  uncertainties  abound  in  the  experimental  procedure.  The 
result  is  that  the  sharp  bifurcation  of  Figure  8.2.  is  much  more  rounded  (as  discussed  later  in 
Figures  8.6  and  8.7),  and  the  critical  buckling  temperature  where  the  first  natural  frequency  is 
zero  is  no  longer  obvious. 

8.2  Methodology  for  State  Boundary  Uncertainty  Quantification 

The  proposed  methodology  quantifies  the  uncertainty  in  a  boundary  between  two  states.  It  is 
assumed  that  the  transition  from  state  i  to  state  j  is  smooth,  and  that  only  states  i  and  j  are  under 
consideration  in  a  given  analysis.  The  methodology  consists  of  three  basic  steps: 

1.  Formulate  a  model  selection  and  state  boundary  uncertainty  quantification  (SBUQ) 
problem 

a.  Identify  possible  operating  states  and  corresponding  models 

b.  Select  appropriate  variable(s)  given  observations  and  models  for  SBUQ. 

2.  Create  Bayesian  network  of  the  problem 
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a.  Define  graphical  structure  corresponding  to  model  selection  problem 

b.  Prepare  training  data  for  BN 

c.  Estimate  BN  distribution  parameters 
3.  Quantify  state  boundary  uncertainty 

a.  Estimate  probability  of  system  being  in  each  operating  state  with  respect  to  a 
particular  measurement  via  inference 

Features  of  the  methodology  are  illustrated  throughout  this  section  using  the  thin  beam 
problem.  In  this  problem,  a  thin  beam  subjected  to  temperature  load  may  be  in  either  a  pre¬ 
buckled  or  post-buckled  configuration.  Due  to  many  sources  of  uncertainty,  discussed  in  Section 
8.3,  it  is  difficult  to  detennine  the  exact  buckling  temperature  of  the  beam,  so  the  distribution  of 
the  buckling  temperature  is  sought. 

8.2.1  Probabilistic  Model  Selection 

Consider  a  problem  domain  over  the  variables  x  G  (x1;  ...,xn}.  x  consists  of  two  subsets,  the 
set  of  r  observable  variables  x  G  {z1, ... ,  zr}  and  the  set  of  of  p  unobservable  (hidden)  variables 
h  G  (bj, ... ,  bp]  such  that  n  —  r  +  p.  For  example,  zi  could  represent  a  measurable  quantity  such 
as  velocity  or  natural  frequency  while  xi  could  represent  an  unobservable  system  variable  such  as 
an  operating  mode,  health  state,  or  even  measurable  quantities  for  which  no  measurements  are 
available.  A  single  realization  of  z;  is  denoted  by  Zi  3  Z and  a  single  realization  of  x;  by  X;.  3  X 
For  the  model  selection  formulation,  a  variable  M3  h,  represents  a  choice  of  system  model 
(operating  mode)  that  cannot  be  directly  measured  or  observed.  As  such,  the  value  of  M,  i.e.,  M  = 
mi  and  M=  m2,  must  be  inferred  by  considering  the  measured  variables  z. 

The  critical  value  of  measurement  variable  zu  denoted  zf,  is  the  boundary  between  two  states 
(models)  M  =  mj  and  M  =  m2  along  dimension  Zj.  For  example,  zf  could,  refer  to  the  stress  at 
which  a  material  switches  from  elastic  to  plastic  behavior.  Estimating  Pr(M  =  rrij  |  Z)  over  a 
range  of  Zj  presumed  to  contain  zf  provides  the  infonnation  to  estimate  the  distribution  of  zf, 
p(zf).  The  estimation  of  Pr(M  =  m;|z)  may  be  viewed  through  the  Bayes  factor  and 
probabilistic  model  selection. 

Given  observations  Z  and  model  choices  M  =  m  1  and  M  =  m2,  the  ratio  of  posterior  odds  for 
competing  models  is 


Pr(M  =  m1\Z')  Pr(Z\M  —  mj  Pr(M  —  m^) 

Pr{M  —  m2\Z )  Pr(Z\M  —  m2)X  Pr(M  —  m2) 

The  first  tenn  on  the  right  hand  side  is  the  Bayes  factor,  and  the  second  tenn  is  the  ratio  of 
prior  odds.  Thus,  the  Bayes  factor  is  the  ratio  of  prior  to  posterior  odds.  When  Pr(M  =  mQ  = 
Pr(M  =  m2),  the  Bayes  factor  is  the  posterior  odds  ratio.  When  the  posterior  odds  ratio  is  equal 
to  1,  both  models  are  equally  likely.  As  B12  increases,  the  evidence  against  M  =  m2  increases, 
providing  a  basis  for  model  selection. 

A  probabilistic  model  selection  problem  can  be  formulated  where  the  model  choice  is 
included  as  a  discrete  random  variable  [93,94],  Suppose  the  problem  is  to  select  from  amongst  K 
possible  model  choices,  each  with  parameter  vector  Oj,j  —  1  Fet  y  be  a  vector  of 

observable  system  responses.  Using  a  Bayesian  model  specification,  the  prior  probability  of  the 
models  is  Uj  =  Pr(M  =  m.j ),  with  YJj=inj  —  1-  The  model  parameters  Oj  are  conditioned  on  M 
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with  prior  p(0y|M  =  my),  and  the  likelihood  function  for  model  j  is  f(z\dj,M  —  my).  The  joint 
posterior  distribution  is 


p(z,  Bj,  M  =  m;)  =  f{z\Oj,  M  =  my) 


K 


f[p(0i\M  = 

i= 1 


my) 


TT; 


(8.5) 


The  beam  buckling  problem  of  Section  8.1  considers  two  possible  configurations  (K  =  2),  i.e., 
pre-buckled  and  post-buckled.  The  models  for  the  first  and  third  natural  frequency  in  each 
configuration  have  their  own  set  of  parameters  and  prior  probabilities,  so  we  have  61  = 
{d{1, 9{3},  d2  =  ni  =  [n 2.1’n2.3}’  anc*  71 1  =  {yrf1, 7rp}.  Where  %  and  n2  depend 

on  the  temperature,  T,  thus  nj  =  p(M  =  my|r).  The  value  of  M,  a  post-buckled  or  pre-buckled 
model,  is  not  always  observable.  The  measured  values  Z  consist  of  the  first  and  third  natural 
frequencies  as  well  as  the  temperature,  T.  Although  M  may  be  observable  in  some  states  of  the 
system  (e.g.  when  buckling  is  clearly  visible  via  the  naked  eye  or  from  strain  gages)  and  included 
in  Z,  this  does  not  necessarily  occur  in  the  region  of  interest. 

The  quantity  of  interest  for  model  selection  is  the  posterior  probability  of  p(M  =  my  |  Oj ,  Z), 


p(M  =  my  Bj, z) 


f{z\Bj,M  -  my)[nr=iP(fli|M  -  my)]?ry 
XLi/(zl0/c<M  =  mkmUv(.Ot\M  =  mk)]nk 


(8.6) 


which  is  also  an  element  in  the  posterior  odds  ratio  of  Eq.  (8.4).  In  terms  of  the  beam  buckling 
problem,  Eq.  (8.6)  gives  the  probability  that  the  beam  is  in  either  the  pre-buckled  or  post-buckled 
configuration.  Given  N  samples  from  the  posterior  distribution,  the  posterior  probability  of  M  = 
ny  is  approximated  by 


„  ,  .  #  samples  where  M  =  m.- 

p(M  =  my  |  z)  = - - - - - L  (8. 7) 

The  formulation  of  Eq.  (8.6)  can  be  expressed  using  a  Bayesian  network  (BN),  as  discussed  in 
the  next  section. 

8.2.2  Bayesian  Network  Formulation  of  Probabilistic  Model  Selection 

A  Bayesian  network  (BN)  represents,  graphically  and  mathematically,  the  joint  probability 
distribution  over  a  set  of  variables,  x  G  (xx  ...xn}.  A  BN  is  a  directed  acyclic  graph  (DAG):  the 
edges  (arrows)  specify  the  dependence  structure  of  x,  which  is  represented  by  nodes.  The  acyclic 
requirement  provides  that  starting  from  node  i,  there  exists  no  path  in  the  network  to  return  to 
node  i.  The  joint  distribution  may  be  written  in  factored  fonn  as 
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Figure  8.3.  Bayesian  network  of  beam  buckling  system 

71 

POl,  -*n)  =  |^[  p(x£  |77|)  (8  8) 

i=l 

where  n£  is  the  set  of  parent  variables  for  x£,  i.e.,  the  variables  on  which  x£  is  immediately 
dependent,  and  p(x£|rij)  is  the  conditional  probability  distribution  (CPD)  for  x£ .  If  n£  =  0, 
p(xj|IIj)  =  p(xj)  is  a  degenerate  conditional  probability  distribution,  but  in  this  paper  the  term 
CPD  shall  include  both  regular  and  degenerate  CPDs. 

The  Bayesian  network  of  the  thin  beam  system  is  shown  in  Figure  8.3  and  the  conditional 
probability  distributions  are  summarized  in  Table  8.1.  The  temperature  load,  T,  is  considered  a 
deterministic  input,  indicated  by  the  square  node.  It  is  assumed  that  T  will  always  be  available 
(this  assumption  could  be  relaxed  by  assuming  a  distribution  for  T).  The  buckling  state,  b,  is  a 
binomial  logistic  regression  for  predicting  whether  the  beam  is  either  pre-  or  post-buckled 
depending  upon  T.  The  first  and  third  natural  frequencies,  f\  and  /?,  are  modeled  as  Gaussian 
processes  based  on  T.  Since  the  natural  frequencies  also  depend  upon  the  buckling  state,  b,  there 
is  a  Gaussian  process  model  for  each  state  for  both  f\  and/?.  Thus,  the  BN  contains  a  total  of  4 
Gaussian  process  models.  Parameterization  of  the  BN  is  discussed  in  the  next  section.  Inference 
in  the  BN  of  Figure  8.3  is  performed  via  likelihood  weighting,  a  form  of  importance  sampling 
that  is  suitable  for  low  dimensional  networks  [28]. 


Table  8.1.  Bayesian  network  variables 


Variable 

Symbol 

CPD 

Cardinality 

Temperature 

T 

Deterministic  Input 

Continuous 

Buckling 

b 

Logistic  Regression 

Binary 

1st  natural  frequency 

fi 

Gaussian  Process  x  2 

Continuous 

3Id  natural  frequency 

fi 

Gaussian  Process  x  2 

Continuous 

8.2.3  Handling  Uncertainty  in  Training  Data  Near  the  Boundary 

Given  that  the  state  boundary  itself  is  uncertain,  data  obtained  near  the  boundary  may  not 
belong  to  an  obvious  state  and  thus  have  no  class  label  associated  with  it.  Such  data  is  called 
unlabeled  data.  This  is  problematic  because  estimating  model  training/calibration  requires 
complete  data  (i.e.,  no  missing  values).  Even  so,  unlabeled  data  may  contain  important 
information  about  system  behavior,  and  has  been  shown  under  appropriate  modeling 
assumptions  to  improve  the  Fisher  information  metric  [96]  (a  measure  of  the  value  of  training 
data)  of  the  probability  model  [97].  For  example,  the  4-tuple  of  data  (T,bJ 1/3)  may  not  have  a 
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value  of  b,  but  still  have  values  for/i  and /  that  can  be  used  for  analysis.  This  situation  occurs 
near  the  state  transition  boundary  because  it  is  difficult  to  pinpoint  the  exact  moment  when 
buckling  occurs  on  a  vibrating  beam  specimen. 


Finding  The  Labeling  Boundary 


Figure  8.4.  The  labeling  boundary  is  a  weighted  average  of  the  temperatures  at  which  the  / 
and /  mean  functions  are  at  their  respective  minimums 

In  the  thin  beam  problem,  the  state  transition  boundary  corresponds  to  the  temperature  at 
which  the  natural  frequencies  are  at  their  minimum.  Figure  8.4  shows  the  mean  Gaussian  Process 
(GP)  model  responses  (with  AT  as  the  independent  variable)  constructed  from  the  training  data 
for/i  and  /.  No  labels  for  pre-  or  post-buckling  are  required  to  fit  these  GP  models.  The 
threshold  AT  value  for  classification  is  unclear  because,  as  seen  in  Figure  8.4,  the  mean  curves 
for/i  and/3  have  minimum  values  at  different  temperatures,  (A T[  and  A//.  This  is  an  indicator 
of  the  uncertainty  in  the  buckling  temperature  and  the  reason  for  being  unable  to  confidently 
label  the  data  in  the  first  place.  To  account  for  the  uncertainty  in  the  training  data,  a  variance 
weighted  labeling  threshold  A Tthresh  —  w1A7’1'  +  w3AT3  is  estimated,  where  wl  —  l  — 
of'  /  (ai'  +  °i')  and  w3  =  1  —  erf'  /  (of'  +  of ') .  the  variances  at  both  minimums  (of'  and 
of')  are  calculated  from  the  GP  models.  The  training  data  can  now  be  labeled  as  either  pre-  or 
post-buckled  using  the  class  boundary,  and  the  4  GP  models  describing  f  and  /  in  the  pre-  and 
post-buckled  states  can  now  be  trained.  These  GP  models  may  be  built  with  standard  covariance 
functions,  such  as  the  Matem  class  and  trained  using  maximum  likelihood  estimation  or 
Bayesian  methods. 

Even  with  an  estimate  of  the  labeling  threshold,  these  GP  models  include  little  information 
about  the  confidence  in  the  labels  assigned  to  the  training  data.  This  labeling  uncertainty  is 
accounted  for  in  the  BN  by  the  logistic  regression:  the  likelihoods  of  the/  and /  predictions  will 
be  combined  with  the  probability  of  being  pre-  or  post-buckled  during  Bayesian  updating  in  the 
BN,  as  determined  by  a  logistic  regression,  represented  by  b  in  the  Bayesian  network  of  Figure 
8.3.  The  logistic  regression  computes  the  probability  of  being  either  pre-  or  post-buckled  as  a 
function  of  temperature.  The  probability  effectively  weights  the  confidence  in  the  responses  from 
the  GP  model  at  a  given  temperature. 
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The  regression  of  a  training  set  of  m  labeled  examples  {(x^.y^), (x(-rn\y(-7n'))},  where 
the  discrete  variable  y(l)  has  k  possible  outcomes  and  x  (a:>  G  Mn+1  with  xo  =  1  as  the  intercept 
may  be  described  by 


■p(y(t)  —  l|xW;  oy 

~ed[xCl)' 

he(xw)  = 

p(yd)  —  2\x^l\  0) 

1 

e0]A(l) 

"aiT 

II 

^7-1 

-e9kx(0 . 

with  parameters  0  =  61,62,  ■■■ ,  0k and  normalizing  term  1  /Ylj= i  e9J x() .  If  k  =  2,  this  is  a  logistic 
regression,  and  if  k  >  2,  the  regression  is  termed  a  multinomial  logistic  regression  or  softmax 
regression  [99]. 

Typically,  the  data  for  training  the  logistic  regression  would  consist  of  labeled  data  where  the 
gradual  change  in  probability  of  observing  the  different  classes  is  reflected  in  the  training  data. 
However,  since  the  training  data  was  labeled  using  an  abrupt  boundary,  a  logistic  regression 
fitted  to  the  training  data  will  reflect  this  by  having  a  sharp  boundary  between  classes.  The 
gradual  transition  between  classes  expected  to  arise  due  to  experimental  uncertainties  will  not  be 
present.  However,  a  process  called  L2  regularization  (or  variously  Tikhonov  [106]  regularization 
or  ridge  regression  [107])  provides  a  means  of  controlling  the  sharpness  of  the  class  boundary.  L2 
regularization  is  often  used  to  limit  model  complexity  by  adding  a  weight  decay  tenn  to  the  cost 
function  being  minimized  during  regression.  Increasing  the  size  of  the  penalty  has  the  effect  of 
smoothing  the  model.  The  weight  decay  tenn  is  often  found  via  cross-validation  to  find  the 
model  that  achieves  the  best  predictive  perfonnance.  Equation  (8.10)  shows  the  cost  function  for 
logistic  regression  with  a  weight  decay  term  with  parameter  X  added  on  the  right  hand  side. 
Equation  (8.10)  is  minimized  for  a  given  X  to  find  the  optimal  parameters  0. 


y(0)  =  — 

m 


m  k 

II 

i= 1  7=1 


/{y(0  =  j^log 


g T 

e  j 


yk  pOfxW 


ill* 

i= 1  ;=0 


(8.10) 


In  this  study,  the  weight  decay  tenn  is  varied  until  the  resulting  logistic  regression  fits  the 
expert’s  beliefs  about  the  probability  of  being  in  the  post-buckling  condition  at  a  particular 
temperature  (e.g.,  1%  probability  at  A  T  =  7°  C).  The  expert  would  make  this  judgment 
qualitatively  or  by  employing  heuristics,  the  goal  being  to  capture  the  region  of  classification 
uncertainty.  Then,  by  finding  the  temperature  corresponding  to  50%  probability,  symmetry  can 
be  used  to  find  the  boundary  on  the  other  side,  as  the  logistic  regression  assumes  the  log  odds  of 
the  probability  varies  linearly.  This  procedure  provides  a  means  to  fuse  expert  opinion  into  the 
methodology. 
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8.2.4  Quantification  of  State  Boundary  Uncertainty 

With  a  trained  BN  available,  the  model  variable,  M,  can  be  inferred.  Given  some  set  of 
measurements,  Z,  that  does  not  include  the  model  variable  M,  the  inference  algorithm  is  run  with 
N  samples  to  estimate  p(M  =  m^Z)  for  all  j.  Of  course,  Z  represents  only  one  instantiation  of 
the  measured  variables.  Thus,  p(M  =  m^Z)  is  only  a  point  estimate  of  probabilities  of  M  —  mj. 
To  determine  the  distribution  of  the  boundary,  p(zf),  the  value  of  p(M  =  m;|z)  must  be 
computed  via  inference  over  a  range  of  values  for  Zj .  While  Zj  values  are  selected  to  represent  a 
range  of  zt,  no  selection  criteria  is  extended  to  z_,;  (-1  i  refers  to  all  indices  except  /).  Figure  8.5 
illustrates  how  the  probability  of  M  being  in  one  of  two  states  changes  with  a  variable  Zj  and  the 
resulting  distribution  of  zf . 


Figure  8.5.  State  boundary  uncertainty 

Left:  Pr(M  =  mi)  and  Pr(M  =  m2)  vary  as  the  boundary  x  =  0  is  approached.  Right:  The 
distribution  of  the  critical  z  value  that  marks  the  boundary  between  M  =  mi  and  M  =  m2. 

In  the  beam  buckling  problem,  p(zf)  is  the  distribution  of  the  buckling  temperature,  which  is 
estimated  by  evaluating  the  probability  of  the  model  being  in  the  pre-buckled  or  post-buckled 
state  at  predetennined  temperature  values.  Temperature  is  the  Zj  variable  and  the  first  and  third 
natural  frequencies,  f\  and  are  the  z_,;  variables,  which  are  allowed  to  vary. 

8.2.5  Methodology  Summary 

A  methodology  was  developed  for  the  quantification  of  state  boundary  uncertainty.  First,  a 
model  selection  problem  is  developed.  The  model  selection  problem  is  then  formulated  as  a  BN. 
Data  near  the  boundary  is  classified  using  an  optimization  algorithm.  The  BN  may  then  be 
trained.  Finally,  inference  is  perfonned  in  the  BN  for  a  range  of  values  of  the  variable  of  interest. 
A  probability  distribution  of  the  state  boundary  is  then  constructed.  Next,  a  case  study  is  used  to 
demonstrate  the  methodology. 

8.3  Case  Study:  Thin  Beam  Buckling 

Section  8.2  outlined  a  methodology  for  quantifying  the  uncertainty  in  the  state  boundary, 
which  is  applied  in  this  section  for  the  pre-  and  post-buckled  states  of  a  thin  beam  subjected  to  a 
thermal  load. 

8.3.1  Experimental  Procedure 

Due  to  the  difficulty  of  observing  buckling  temperature  directly  with  any  reasonable 
confidence  (particularly  in  the  region  of  interest),  Perez  et  al.  developed  an  experiment  to 
measure  the  natural  frequencies  of  a  thin  steel  beam. 
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The  thin  steel  beam  ( E  =  200  GPa,  p  =  7860  kg/m  )  was  mounted  on  a  clamped-clamped 
support  with  two  bolts  tightened  with  a  torque  wrench  on  either  end.  The  beam  dimensions  were 
22.86  cm  (x-direction)  x  1.27  cm  (y-direction)  x  0.0788  cm  (z-direction).  The  straightest  beam 
specimen  available  was  selected  to  minimize  the  initial  curvature  of  the  beam.  The  beam  was 
subjected  to  vibration,  Fz,  in  the  z-axis  from  a  magnetic  driver,  which  was  supplied  a  random 
noise  signal.  The  driver  was  positioned  near  the  support  to  minimize  the  oscillation  amplitude. 
The  beam  deflection  was  measured  with  a  laser  vibrometer  at  the  center  of  the  beam.  From  these 
measurements,  the  transient  power  spectral  density  (PSD)  was  computed.  Time  histories  of 
natural  frequencies  were  then  extracted  from  the  PSDs. 

Thennocouples  were  attached  to  the  beam  to  monitor  the  temperature.  The  beam  was  heated 
to  a  post-buckled  state  and  held  until  steady  state  temperature  was  achieved,  then  allowed  to 
cool.  As  the  beam  cooled,  natural  frequencies  and  temperatures  were  recorded  at  specific 
intervals.  The  heating  and  cooling  process  was  repeated  five  times  with  the  same  clamped  beam 
specimen  to  minimize  variability  in  boundary  conditions  between  experiments. 

8.3.2  Uncertainty  in  Buckling  Temperature  Identification 

Many  uncertainties  arise  out  of  the  experimental  process  through  the  combination  of 
laboratory  environment  uncertainties,  test  equipment  limitations,  and  computational  errors. 
These  uncertainties  result  in  behavior  that  deviates  from  the  predictions  of  idealized  models. 
Deflection  measurements  from  the  laser  vibrometer  may  be  erroneous  if  the  laser  is  not  reflected 
directly  back  at  the  source,  which  contains  a  sensor  for  measuring  distance.  This  may  happen  if 
the  laser  is  not  pointed  directly  at  an  inflection  point  and  the  deflection  is  large. 

Temperature  readings  contain  uncertainties  due  to  changing  ambient  air  conditions  (largely 
due  to  climate  control  systems  in  the  laboratory)  and  heating  of  the  beam,  which  is  difficult  to 
control.  The  ambient  temperature  was  between  13  and  14  degrees  Celsius  for  each  of  the  5  tests. 
Uncalibrated  or  damaged  thennocouples  can  produce  erroneous  data.  Furthermore,  due  to 
equipment  limitations,  only  initial  and  final  temperature  values  are  measured  at  3  locations  on 
the  beam  (left,  middle,  and  right).  A  quadratic  cooling  model  is  used  to  calculate  intermediate 
temperatures.  Zero  mean  Gaussian  noise  with  a  coefficient  of  variation  of  0.01  is  added  to  each 
initial  temperature  to  simulate  the  effect  of  measurement  noise  before  fitting  the  cooling  model. 
Synchronization  with  deflection  (and  subsequently  frequency)  measurements  is  difficult  to 
achieve,  as  the  acquisition  systems  for  each  measurement  are  independently  operated.  Noise  is 
added  back  to  the  temperature  measurements  by  adding  zero  mean  Gaussian  noise  with  a 
standard  deviation  of  0.025,  a  value  experimentally  detennined  to  approximate  the  noise  in 
temperature  measurements. 

The  first  and  third  natural  frequencies  (f\  and  f$)  are  extracted  from  the  power  spectral 
densities  (PSDs)  of  the  beam  response.  However,  the  PSDs  are  only  approximations  that  are 
based  on  a  window  of  transient  data.  The  time  length  of  the  window  affects  the  smoothness  of 
the  PSD.  A  peak  detection  algorithm  is  used  to  find  the  natural  frequencies  within  the  PSD. 

Other  sources  of  uncertainty  (e.g.,  initial  stress  state  and  boundary  conditions  of  the  beam, 
initial  curvature  and  geometry  of  the  beam,  locations  of  heat  lamps)  make  repeatability  of  the 
experiment  difficult. 

8.3.3  Experimental  Data 

Figure  8.6  shows  experimental  data  for  the  first  four  iterations  of  the  experiment.  Each 
iteration  consists  of  heating  and  cooling  the  beam.  The  same  beam  and  boundary  conditions  are 
used  for  each  iteration.  The  sharp  bifurcation  seen  in  the  idealized  frequency  response  in  Figure 
8.2  is  blunted  in  both  Figures  8.6  and  8.7.  It  was  not  possible  to  experimentally  obtain  0  Hz 
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natural  frequencies  due  to  uncertainties.  It  is  difficult  to  ascertain  whether  the  beam  is  in  a  pre-  or 
post-buckled  condition  between  about  3.5  and  6°  C  above  ambient  temperature  because  the 
results  do  not  converge  on  a  clear  buckling  temperature.  Furthermore,  the  first  natural  frequency 
suggests  a  lower  buckling  temperature  than  the  third  natural  frequency.  Thus,  it  is  difficult  to 
confidently  label  data  in  this  region.  The  classification  boundary  determined  by  finding  the 
minimum  frequencies  via  optimization  is  shown  in  Figure  8.4. 
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Figure  8.6.  Four  experiments  with  the  expert’s  99%  confidence  bounds  for  data  labeling 
8.3.4  Bayesian  Network  and  SBUQ 

The  BN  of  the  system  is  shown  in  Figure  8.3  with  parameters  described  in  Table  8.1.  The  GP 
models  are  built  using  the  Gaussian  Processes  for  Machine  Learning  (GPML)  toolbox  for  Matlab 
and  have  Matern  covariance  functions  [108].  Figure  8.7  shows  GP  model  fits  with  +/-  3er 
(standard  deviation)  bounds  using  the  first  four  sets  of  experimental  data.  The  data  has  been 
labeled  using  the  procedure  discussed  in  Section  8.2,  resulting  in  a  different  set  of  training  data 
for  each  of  the  four  GP  models.  It  can  be  seen  that  each  model  has  regions  of  lower  variance, 
where  it  has  training  data  (training  data  is  color  coded  for  each  model),  and  regions  of  larger 
variance,  where  the  model  extrapolates.  Similar  variances  for  the  models  around  the  buckling 
temperature  indicate  that  neither  model  is  strongly  preferred  in  this  region. 

To  fit  the  logistic  regression  of  the  discrete  variable  b  on  temperature,  T,  the  analyst  selects  a 
temperature  value  where  they  have  99%  confidence  that  the  beam  is  pre-buckled,  based  on  the 
first  four  tests.  This  value  is  3.25°  C  and  is  indicated  by  the  left  thick  black  line  in  Figure  8.8.  An 
upper  value  of  the  region  of  uncertainty  is  determined  via  symmetry  as  7.33°  C.  The  penalty 
parameter,  2,  is  detennined  via  optimization.  The  teal  line  in  Figure  8.8  shows  the  resulting 
logistic  regression  from  the  first  four  tests.  The  more  uncertainty  in  the  data,  the  greater  the 
distance  between  thick  black  lines. 
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GP  Model  Fits  and  Training  Data  with  +/-3  c  bounds 
Original  Signal 


Figure  8.7.  GP  model  fits,  training  data,  and  3<r  prediction  bounds. 
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Figure  8.8.  Normal  CDFs  of  original  data  fitted  via  least  squares.  Data  points  are  shown. 
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Figure  8.9.  PDFs  of  buckling  temperature  corresponding  to  fitted  CDFs. 

The  teal  line  in  Figure  8.8  corresponds  to  the  post-buckling  probability  calculated  using  the 
logistic  regression  with  temperature  values  from  test  5.  The  colors  blue,  green,  and  red  refer  to 
scenarios  where  in  addition  to  T,  measurements  of and  are  included  in  the  analysis, 

respectively.  These  estimated  probabilites  of  post-buckling,  denoted  by  an  ‘x’  in  Figure  8.8,  are 
calculated  as  in  Section  8.2. 

The  estimated  probabilities  in  Figure  8.8  do  not  necessarily  result  in  an  increasing  function 
due  to  noise  in  the  training  data  and  errors  in  estimating  the  GP  parameters.  At  a  given  value  of 
T,  the  probability  of  post-buckling  is  a  distribution,  and  test  5  only  provides  point  estimates  of f\ 
and  fi,  resulting  in  only  one  sample  of  the  post-buckling  probability.  Because  there  is  no 
guarantee  that  this  probability  increases  with  T,  the  probability  of  post-buckling  as  a  function  of 
temperature  cannot  be  interpreted  as  a  CDF,  as  would  be  desired.  However,  in  situations  where 
decreasing  probabilites  can  be  attributed  to  noise  or  outliers  in  the  data  that  can  be  ignored,  a 
couple  options  are  available  to  obtain  a  CDF.  1)  A  parametric  CDF  may  be  fit  to  the  probabilities 
or  2)  the  probabilities  may  be  smoothed/adjusted  to  ensure  that  the  probabilites  are  always 
increasing  to  construct  an  empirical  CDF.  The  smooth  lines  in  Figure  8.8  correspond  to  least 
squares  estimates  of  Gaussian  CDFs  to  the  post-buckling  probabilities.  The  curves  appear  to  fit 
the  data  quite  well,  qualitatively.  Quantitative  goodness-of-fit  tests  are  problematic  because  they 
assume  a  true  empirical  CDF  is  being  tested  or  can  be  skewed  by  a  single  outlier.  Figure  8.9 
shows  the  Gaussian  PDFs  corresponding  to  the  CDFs  of  Figure  8.8.  The  Gaussian  distribution 
parameters  are  shown  in  Table  8.2.  The  standard  deviation  and  coefficient  of  variation  of  the 
buckling  temperature  decrease  as  more  measurement  variables  are  considered,  indicating  an 
increase  in  precision  with  information.  It  is  not  possible  to  assess  the  solution  bias,  as  the  true 
buckling  temperature  is  unknown. 
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Table  8.2.  Gaussian  distribution  parameters. 


Scenario 

Mean  (°C) 

Standard  Deviation 

Coefficient  of 
Variation 

T,  fi,  and  f3  measured 

5.07 

0.08 

0.016 

T  and  /)  measured 

5.13 

0.21 

0.041 

T  and  /i  measured 

5.10 

0.49 

0.096 

T  measured 

5.28 

0.86 

0.16 

8.3.5  Results  and  Discussion 

Experimental  and  computational  uncertainties  impede  confident  labeling  of  the  system  state. 
This  results  in  a  semi-supervised  learning  problem,  where  some  of  the  training  examples  can  be 
classified,  while  others,  in  the  region  of  interest  near  the  boundary,  cannot.  The  proposed 
methodology  fuses  expert  opinion  with  both  discriminative  and  generative  classification  methods 
in  a  Bayesian  network  framework  to  estimate  uncertainty  in  the  critical  buckling  temperature  in 
the  presence  of  such  classification  uncertainties.  Discriminative  methods  (e.g.,  logistic 
regression)  make  weaker  assumptions  about  underlying  distributions  but  require  more  training 
data  than  generative  methods  (e.g.,  Gaussian  discriminant  analysis),  which  trade  data  for  stronger 
assumptions  about  underlying  distributions. 

In  the  example  problem,  the  logistic  regression  is  the  discriminative  classification  method  and 
the  GP  models  represent  the  generative  method.  The  weight  decay  parameter  A  that  results  from 
regularization  of  the  logistic  regression  is  a  convenient  way  to  account  for  expert  opinion.  In  the 
beam  example,  A  was  chosen  such  that  the  logistic  regression  matched  where  the  analyst  was  at 
least  99%  in  the  labels  of  the  experimental  data.  Without  this  regularization,  the  logistic 
regression  would  simply  show  a  sharp  change  in  probability  at  the  temperature  Top,  which  was 
detennined  by  using  the  temperature  where  the  frequencies  are  at  a  minimum  for  labeling  the 
data.  This  would  result  in  overconfidence  in  the  estimate  of  the  buckling  state. 

Uncertainty  in  classification  based  on  temperature  passes  from  the  regularized  logistic 
regression  to  the  GP  models  during  Monte  Carlo  inference.  Uncertainties  arise  in  the  collection 
and  estimation  of  temperature  data  and  in  computational  errors  in  the  frequency  data  due  to 
approximating  the  FFT  from  velocity  measurements  collected  over  a  window  of  time.  These 
uncertainties  manifest  themselves  as  variance  in  GP  predictions,  which  is  then  reflected  in 
decreased  likelihoods  for  frequency  measurements.  Further,  data  points  near  the  buckling 
boundary  may  have  initially  been  assigned  to  the  wrong  GP  model,  causing  the  wrong  GP  model 
to  produce  larger  likelihoods  near  the  boundary.  Fikelihoods  resulting  from  extrapolation  outside 
the  training  data  may  have  limited  value.  However,  the  logistic  regression  helps  guard  against 
these  issues  by  producing  varying  the  amount  of  samples  associated  with  each  buckling  state 
according  to  expert  opinion.  Ultimately,  the  estimated  bias  and  variance  of  the  critical  buckling 
temperature  distribution  both  depend  on  the  quantity  and  the  quality  of  the  experimental  data, 
modeling  assumptions,  and  computational  procedures. 

The  example  problem  showed  how  uncertainty  in  the  buckling  temperature  of  a  beam  may  be 
quantified,  and  can  result  in  a  parametric  distribution  for  the  buckling  temperature.  The  reduction 
in  uncertainty  by  providing  additional  measurement  data  to  the  analysis  is  evident  from  the 
decreasing  variance  of  the  distributions.  Reduction  of  uncertainty  is  not  guaranteed,  however.  If 
the  buckling  temperature  were  insensitive  to  a  particular  measurement,  there  would  be  little 
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effect  on  the  buckling  temperature  distribution.  Or,  a  faulty  sensor  could  increase  uncertainty  in 
the  buckling  temperature. 

It  should  be  noted  that  the  methodology  was  demonstrated  using  particular  statistical 
techniques  (logistic  regression,  Gaussian  process  models,  likelihood  weighting  inference). 
However,  these  need  not  be  the  only  options.  Many  other  classification  techniques  could  take  the 
place  of  logistic  regression,  e.g.,  support  vector  machines,  neural  networks,  Gaussian  processes. 
The  GP  models  for  the  natural  frequencies  could  instead  be  relevance  vector  machines  or 
polynomial  chaos  models  or  some  other  response  surface  method.  Likelihood  weighting 
inference  could  be  replace  with  a  Gibbs  sampling  scheme. 

8.4  Summary  of  State  Boundary  Uncertainty  Quantification 

The  ability  to  assess  uncertainty  in  state  boundaries  is  of  particular  importance  for  complex 
systems  that  have  many  discrete  model  choices.  This  section  discussed  a  methodology  for 
quantifying  the  uncertainty  in  a  boundary  between  two  system  states,  namely  buckling  in  a  thin 
beam.  First,  the  system  is  cast  in  a  Bayesian  network  framework  containing  a  logistic  regression 
classifier  and  Gaussian  process  surrogate  models.  The  BN  is  trained  using  data  labeled  via  an 
optimization  procedure.  Inference  in  the  Bayesian  network  is  performed  using  likelihood 
weighting.  Uncertainty  in  the  correct  model  of  the  system  is  quantified  over  a  range  of  values  of 
a  measured  variable  to  obtain  the  distribution  of  the  state  boundary.  A  Gaussian  CDF  is  fit  to  the 
probability  of  the  beam  being  post-buckled  as  a  function  of  temperature,  and  it  is  seen  that  the 
methodology  provides  more  precise  estimates  as  more  measurements  are  made  available.  The 
methodology  was  demonstrated  using  experimental  data  from  a  thin  beam  experiment. 
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